Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

colvaratoms.C

Go to the documentation of this file.
00001 #include "colvarmodule.h"
00002 #include "colvarparse.h"
00003 #include "colvaratoms.h"
00004 
00005 
00006 // member functions of the "atom" class depend tightly on the MD interface, and are
00007 // thus defined in colvarproxy_xxx.cpp
00008 
00009 // in this file only atom_group functions are defined
00010 
00011 
00012 // Note: "conf" is the configuration of the cvc who is using this atom group;
00013 // "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
00014 cvm::atom_group::atom_group (std::string const &conf,
00015                              char const        *key)
00016   : b_center (false), b_rotate (false), b_user_defined_fit (false),
00017     ref_pos_group (NULL),
00018     b_fit_gradients (false),
00019     noforce (false)
00020 {
00021   cvm::log ("Defining atom group \""+
00022             std::string (key)+"\".\n");
00023   // real work is done by parse
00024   parse (conf, key);
00025 }
00026 
00027 
00028 cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
00029   : b_dummy (false), b_center (false), b_rotate (false), 
00030     ref_pos_group (NULL), b_fit_gradients (false), 
00031     noforce (false)
00032 {
00033   this->reserve (atoms.size());
00034   for (size_t i = 0; i < atoms.size(); i++) {
00035     this->push_back (atoms[i]);
00036   }
00037   total_mass = 0.0;
00038   for (cvm::atom_iter ai = this->begin();
00039        ai != this->end(); ai++) {
00040     total_mass += ai->mass;
00041   }
00042 }
00043 
00044 
00045 cvm::atom_group::atom_group()
00046   : b_dummy (false), b_center (false), b_rotate (false), 
00047     ref_pos_group (NULL), b_fit_gradients (false), 
00048     noforce (false)
00049 {
00050   total_mass = 0.0;
00051 }
00052 
00053 
00054 cvm::atom_group::~atom_group()
00055 {
00056   if (ref_pos_group) {
00057     delete ref_pos_group;
00058     ref_pos_group = NULL;
00059   }
00060 }
00061 
00062 
00063 void cvm::atom_group::add_atom (cvm::atom const &a)
00064 {
00065   if (b_dummy) {
00066     cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n");
00067   } else {
00068     this->push_back (a);
00069     total_mass += a.mass;
00070   }
00071 }
00072 
00073 
00074 void cvm::atom_group::parse (std::string const &conf,
00075                              char const        *key)
00076 {
00077   std::string group_conf;
00078 
00079   // save_delimiters is set to false for this call, because "conf" is
00080   // not the config string of this group, but of its parent object
00081   // (which has already taken care of the delimiters)
00082   save_delimiters = false;
00083   key_lookup (conf, key, group_conf, dummy_pos);
00084   // restoring the normal value, because we do want keywords checked
00085   // inside "group_conf"
00086   save_delimiters = true;
00087 
00088   if (group_conf.size() == 0) {
00089     cvm::fatal_error ("Error: atom group \""+
00090                       std::string (key)+"\" is set, but "
00091                       "has no definition.\n");
00092   }
00093 
00094   cvm::increase_depth();
00095 
00096   cvm::log ("Initializing atom group \""+std::string (key)+"\".\n");
00097 
00098   // whether or not to include messages in the log
00099   // colvarparse::Parse_Mode mode = parse_silent;
00100   // {
00101   //   bool b_verbose;
00102   //   get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
00103   //   if (b_verbose) mode = parse_normal;
00104   // }
00105   colvarparse::Parse_Mode mode = parse_normal;
00106 
00107   {
00108     //    std::vector<int> atom_indexes;
00109     std::string numbers_conf = "";
00110     size_t pos = 0;
00111     std::vector<int> atom_indexes;
00112     while (key_lookup (group_conf, "atomNumbers", numbers_conf, pos)) {
00113       if (numbers_conf.size()) {
00114         std::istringstream is (numbers_conf);
00115         int ia;
00116         while (is >> ia) {
00117           atom_indexes.push_back (ia);
00118         }
00119       }
00120 
00121       if (atom_indexes.size()) {
00122         this->reserve (this->size()+atom_indexes.size());
00123         for (size_t i = 0; i < atom_indexes.size(); i++) {
00124           this->push_back (cvm::atom (atom_indexes[i]));
00125         }
00126       } else
00127         cvm::fatal_error ("Error: no numbers provided for \""
00128                           "atomNumbers\".\n");
00129 
00130       atom_indexes.clear();
00131     }
00132   }
00133 
00134   {
00135     std::string range_conf = "";
00136     size_t pos = 0;
00137     while (key_lookup (group_conf, "atomNumbersRange",
00138                        range_conf, pos)) {
00139 
00140       if (range_conf.size()) {
00141         std::istringstream is (range_conf);
00142         int initial, final;
00143         char dash;
00144         if ( (is >> initial) && (initial > 0) &&
00145              (is >> dash) && (dash == '-') &&
00146              (is >> final) && (final > 0) ) {
00147           for (int anum = initial; anum <= final; anum++) {
00148             this->push_back (cvm::atom (anum));
00149           }
00150           range_conf = "";
00151           continue;
00152         }
00153       }
00154 
00155       cvm::fatal_error ("Error: no valid definition for \""
00156                         "atomNumbersRange\", \""+
00157                         range_conf+"\".\n");
00158     }
00159   }
00160 
00161   std::vector<std::string> psf_segids;
00162   get_keyval (group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode);
00163   for (std::vector<std::string>::iterator psii = psf_segids.begin();
00164        psii < psf_segids.end(); psii++) {
00165 
00166     if ( (psii->size() == 0) || (psii->size() > 4) ) {
00167       cvm::fatal_error ("Error: invalid segmend identifier provided, \""+
00168                         (*psii)+"\".\n");
00169     }
00170   }
00171 
00172   {
00173     std::string range_conf = "";
00174     size_t pos = 0;
00175     size_t range_count = 0;
00176     std::vector<std::string>::iterator psii = psf_segids.begin();
00177     while (key_lookup (group_conf, "atomNameResidueRange",
00178                        range_conf, pos)) {
00179       range_count++;
00180 
00181       if (range_count > psf_segids.size()) {
00182         cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than "
00183                           "values of \"psfSegID\".\n");
00184       }
00185 
00186       std::string const &psf_segid = psf_segids.size() ? *psii : std::string ("");
00187 
00188       if (range_conf.size()) {
00189 
00190         std::istringstream is (range_conf);
00191         std::string atom_name;
00192         int initial, final;
00193         char dash;
00194         if ( (is >> atom_name) && (atom_name.size())  &&
00195              (is >> initial) && (initial > 0) &&
00196              (is >> dash) && (dash == '-') &&
00197              (is >> final) && (final > 0) ) {
00198           for (int resid = initial; resid <= final; resid++) {
00199             this->push_back (cvm::atom (resid, atom_name, psf_segid));
00200           }
00201           range_conf = "";
00202         } else {
00203           cvm::fatal_error ("Error: cannot parse definition for \""
00204                             "atomNameResidueRange\", \""+
00205                             range_conf+"\".\n");
00206         }
00207 
00208       } else {
00209         cvm::fatal_error ("Error: atomNameResidueRange with empty definition.\n");
00210       }
00211 
00212       if (psf_segid.size())
00213         psii++;
00214     }
00215   }
00216 
00217   {
00218     // read the atoms from a file
00219     std::string atoms_file_name;
00220     if (get_keyval (group_conf, "atomsFile", atoms_file_name, std::string (""), mode)) {
00221 
00222       std::string atoms_col;
00223       if (!get_keyval (group_conf, "atomsCol", atoms_col, std::string (""), mode)) {
00224         cvm::fatal_error ("Error: parameter atomsCol is required if atomsFile is set.\n");
00225       }
00226 
00227       double atoms_col_value;
00228       bool const atoms_col_value_defined = get_keyval (group_conf, "atomsColValue", atoms_col_value, 0.0, mode);
00229       if (atoms_col_value_defined && (!atoms_col_value))
00230         cvm::fatal_error ("Error: atomsColValue, "
00231                           "if provided, must be non-zero.\n");
00232 
00233       cvm::load_atoms (atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
00234     }
00235   }
00236 
00237   for (std::vector<cvm::atom>::iterator a1 = this->begin(); 
00238        a1 != this->end(); a1++) {
00239     std::vector<cvm::atom>::iterator a2 = a1;
00240     ++a2;
00241     for ( ; a2 != this->end(); a2++) {
00242       if (a1->id == a2->id) {
00243         if (cvm::debug())
00244           cvm::log ("Discarding doubly counted atom with number "+
00245                     cvm::to_str (a1->id+1)+".\n");
00246         a2 = this->erase (a2);
00247         if (a2 == this->end())
00248           break;
00249       }
00250     }
00251   }
00252 
00253   if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
00254     b_dummy = true;
00255     this->total_mass = 1.0;
00256   } else 
00257     b_dummy = false;
00258 
00259   if (b_dummy && (this->size())) 
00260     cvm::fatal_error ("Error: cannot set up group \""+
00261                       std::string (key)+"\" as a dummy atom "
00262                       "and provide it with atom definitions.\n");
00263 
00264 #if (! defined (COLVARS_STANDALONE))
00265   if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) {
00266     cvm::fatal_error ("Error: no atoms defined for atom group \""+
00267                       std::string (key)+"\".\n");
00268   }
00269 #endif
00270 
00271   if (!b_dummy) {
00272     this->total_mass = 0.0;
00273     for (cvm::atom_iter ai = this->begin();
00274          ai != this->end(); ai++) {
00275       this->total_mass += ai->mass;
00276     }
00277   }
00278 
00279   if (!b_dummy) {
00280     bool enable_forces = true;
00281     // disableForces is deprecated
00282     if (get_keyval (group_conf, "enableForces", enable_forces, true, mode)) {
00283       noforce = !enable_forces;
00284     } else {
00285       get_keyval (group_conf, "disableForces", noforce, false, mode);
00286     }
00287   }
00288 
00289   // FITTING OPTIONS
00290 
00291   bool b_defined_center = get_keyval (group_conf, "centerReference", b_center, false, mode);
00292   bool b_defined_rotate = get_keyval (group_conf, "rotateReference", b_rotate, false, mode);
00293   // is the user setting explicit options?
00294   b_user_defined_fit = b_defined_center || b_defined_rotate;
00295 
00296   get_keyval (group_conf, "enableFitGradients", b_fit_gradients, true, mode);
00297 
00298   if (b_center || b_rotate) {
00299 
00300     if (b_dummy)
00301       cvm::fatal_error ("Error: centerReference or rotateReference "
00302                         "cannot be defined for a dummy atom.\n");
00303 
00304     if (key_lookup (group_conf, "refPositionsGroup")) {
00305       // instead of this group, define another group to compute the fit
00306       if (ref_pos_group) {
00307         cvm::fatal_error ("Error: the atom group \""+
00308                           std::string (key)+"\" has already a reference group "
00309                           "for the rototranslational fit, which was communicated by the "
00310                           "colvar component.  You should not use refPositionsGroup "
00311                           "in this case.\n");
00312       }
00313       cvm::log ("Within atom group \""+std::string (key)+"\":\n");
00314       ref_pos_group = new atom_group (group_conf, "refPositionsGroup");
00315     }
00316 
00317     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
00318 
00319     get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode);
00320 
00321     std::string ref_pos_file;
00322     if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) {
00323 
00324       if (ref_pos.size()) {
00325         cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
00326                           "\"refPositions\" at the same time.\n");
00327       }
00328 
00329       std::string ref_pos_col;
00330       double ref_pos_col_value;
00331       
00332       if (get_keyval (group_conf, "refPositionsCol", ref_pos_col, std::string (""), mode)) {
00333         // if provided, use PDB column to select coordinates
00334         bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);
00335         if (found && !ref_pos_col_value)
00336           cvm::fatal_error ("Error: refPositionsColValue, "
00337                             "if provided, must be non-zero.\n");
00338       } else {
00339         // if not, rely on existing atom indices for the group
00340         group_for_fit->create_sorted_ids();
00341         ref_pos.resize (group_for_fit->size());
00342       }
00343 
00344       cvm::load_coords (ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids,
00345                         ref_pos_col, ref_pos_col_value);
00346     }
00347 
00348     if (ref_pos.size()) {
00349 
00350       if (b_rotate) {
00351         if (ref_pos.size() != group_for_fit->size())
00352           cvm::fatal_error ("Error: the number of reference positions provided ("+
00353                             cvm::to_str (ref_pos.size())+
00354                             ") does not match the number of atoms within \""+
00355                             std::string (key)+
00356                             "\" ("+cvm::to_str (group_for_fit->size())+
00357                             "): to perform a rotational fit, "+
00358                             "these numbers should be equal.\n");
00359       }
00360 
00361       // save the center of geometry of ref_pos and subtract it
00362       center_ref_pos();
00363 
00364     } else {
00365 #if (! defined (COLVARS_STANDALONE))
00366       cvm::fatal_error ("Error: no reference positions provided.\n");
00367 #endif
00368     }
00369 
00370     if (b_fit_gradients) {
00371       group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::atom_pos (0.0, 0.0, 0.0));
00372       rot.request_group1_gradients (group_for_fit->size());
00373     }
00374 
00375     if (b_rotate && !noforce) {
00376       cvm::log ("Warning: atom group \""+std::string (key)+
00377                 "\" will be aligned to a fixed orientation given by the reference positions provided.  "
00378                 "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
00379                 "to its radius of gyration), the optimal rotation and its gradients may become discontinuous.  "
00380                 "If that happens, use refPositionsGroup (or a different definition for it if already defined) "
00381                 "to align the coordinates.\n");
00382       // initialize rot member data
00383       rot.request_group1_gradients (this->size());
00384     }
00385 
00386   }
00387 
00388   if (cvm::debug())
00389     cvm::log ("Done initializing atom group with name \""+
00390               std::string (key)+"\".\n");
00391 
00392   this->check_keywords (group_conf, key);
00393 
00394   cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
00395             cvm::to_str (this->size())+" atoms initialized: total mass = "+
00396             cvm::to_str (this->total_mass)+".\n");
00397 
00398   cvm::decrease_depth();
00399 }
00400 
00401 
00402 void cvm::atom_group::create_sorted_ids (void)
00403 {
00404   // Only do the work if the vector is not yet populated
00405   if (sorted_ids.size())
00406     return;
00407 
00408   std::list<int> temp_id_list;
00409   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
00410     temp_id_list.push_back (ai->id);
00411   }
00412   temp_id_list.sort();
00413   temp_id_list.unique();
00414   if (temp_id_list.size() != this->size()) {
00415     cvm::fatal_error ("Error: duplicate atom IDs in atom group? (found " +
00416                       cvm::to_str(temp_id_list.size()) +
00417                       " unique atom IDs instead of" +
00418                       cvm::to_str(this->size()) + ").\n");
00419   }
00420   sorted_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
00421   return;
00422 }
00423 
00424 void cvm::atom_group::center_ref_pos()
00425 {
00426   // save the center of geometry of ref_pos and then subtract it from
00427   // them; in this way it will be possible to use ref_pos also for
00428   // the rotational fit
00429   // This is called either by atom_group::parse or by CVCs that set
00430   // reference positions (eg. RMSD, eigenvector)
00431   ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0);
00432   std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
00433   for ( ; pi != ref_pos.end(); pi++) {
00434     ref_pos_cog += *pi;
00435   }
00436   ref_pos_cog /= (cvm::real) ref_pos.size();
00437   for (std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
00438        pi != ref_pos.end(); pi++) {
00439     (*pi) -= ref_pos_cog;
00440   }
00441 }
00442 
00443 void cvm::atom_group::read_positions()
00444 {
00445   if (b_dummy) return;
00446 
00447   for (cvm::atom_iter ai = this->begin();
00448        ai != this->end(); ai++) {
00449     ai->read_position();
00450   }
00451 
00452   if (ref_pos_group)
00453     ref_pos_group->read_positions();
00454 }
00455 
00456 void cvm::atom_group::calc_apply_roto_translation()
00457 {
00458   atom_group *fit_group = ref_pos_group ? ref_pos_group : this;
00459 
00460   if (b_center) {
00461     // center on the origin first
00462     cvm::atom_pos const cog = fit_group->center_of_geometry();
00463     for (cvm::atom_iter ai = this->begin();
00464          ai != this->end(); ai++) {
00465       ai->pos -= cog;
00466     }
00467   }
00468 
00469   if (b_rotate) {
00470     // rotate the group (around the center of geometry if b_center is
00471     // true, around the origin otherwise)
00472     rot.calc_optimal_rotation (fit_group->positions(), ref_pos);
00473 
00474     for (cvm::atom_iter ai = this->begin();
00475          ai != this->end(); ai++) {
00476       ai->pos = rot.rotate (ai->pos);
00477     }
00478   }
00479 
00480   if (b_center) {
00481     // align with the center of geometry of ref_pos
00482     for (cvm::atom_iter ai = this->begin();
00483          ai != this->end(); ai++) {
00484       ai->pos += ref_pos_cog;
00485     }
00486   }
00487 }
00488 
00489 void cvm::atom_group::apply_translation (cvm::rvector const &t) 
00490 {
00491   if (b_dummy) return;
00492 
00493   for (cvm::atom_iter ai = this->begin();
00494        ai != this->end(); ai++) {
00495     ai->pos += t;
00496   }
00497 }
00498 
00499 void cvm::atom_group::apply_rotation (cvm::rotation const &rot) 
00500 {
00501   if (b_dummy) return;
00502 
00503   for (cvm::atom_iter ai = this->begin();
00504        ai != this->end(); ai++) {
00505     ai->pos = rot.rotate (ai->pos);
00506   }
00507 }
00508 
00509 
00510 void cvm::atom_group::read_velocities()
00511 {
00512   if (b_dummy) return;
00513 
00514   if (b_rotate) {
00515 
00516     for (cvm::atom_iter ai = this->begin();
00517          ai != this->end(); ai++) {
00518       ai->read_velocity();
00519       ai->vel = rot.rotate (ai->vel);
00520     }
00521 
00522   } else {
00523 
00524     for (cvm::atom_iter ai = this->begin();
00525          ai != this->end(); ai++) {
00526       ai->read_velocity();
00527     }
00528   }
00529 }
00530 
00531 void cvm::atom_group::read_system_forces()
00532 {
00533   if (b_dummy) return;
00534 
00535   if (b_rotate) {
00536 
00537     for (cvm::atom_iter ai = this->begin();
00538          ai != this->end(); ai++) {
00539       ai->read_system_force();
00540       ai->system_force = rot.rotate (ai->system_force);
00541     }
00542 
00543   } else {
00544 
00545     for (cvm::atom_iter ai = this->begin();
00546          ai != this->end(); ai++) {
00547       ai->read_system_force();
00548     }
00549   }
00550 }
00551 
00552 cvm::atom_pos cvm::atom_group::center_of_geometry() const
00553 {
00554   if (b_dummy)
00555     return dummy_atom_pos;
00556 
00557   cvm::atom_pos cog (0.0, 0.0, 0.0);
00558   for (cvm::atom_const_iter ai = this->begin();
00559        ai != this->end(); ai++) {
00560     cog += ai->pos;
00561   }
00562   cog /= this->size();
00563   return cog;
00564 }
00565 
00566 cvm::atom_pos cvm::atom_group::center_of_mass() const
00567 {
00568   if (b_dummy)
00569     return dummy_atom_pos;
00570 
00571   cvm::atom_pos com (0.0, 0.0, 0.0);
00572   for (cvm::atom_const_iter ai = this->begin();
00573        ai != this->end(); ai++) {
00574     com += ai->mass * ai->pos;
00575   }
00576   com /= this->total_mass;
00577   return com;
00578 }
00579 
00580 
00581 void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad)
00582 {
00583   if (b_dummy) return;
00584 
00585   for (cvm::atom_iter ai = this->begin();
00586        ai != this->end(); ai++) {
00587     ai->grad = (ai->mass/this->total_mass) * grad;
00588   }
00589 }
00590 
00591 
00592 void cvm::atom_group::calc_fit_gradients()
00593 {
00594   if (b_dummy) return;
00595 
00596   if ((!b_center) && (!b_rotate)) return; // no fit
00597 
00598   if (cvm::debug())
00599     cvm::log ("Calculating fit gradients.\n");
00600 
00601   atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
00602   group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::rvector (0.0, 0.0, 0.0));
00603 
00604   if (b_center) {
00605     // add the center of geometry contribution to the gradients
00606     for (size_t i = 0; i < this->size(); i++) {
00607       // need to bring the gradients in original frame first
00608       cvm::rvector const atom_grad = b_rotate ?
00609         (rot.inverse()).rotate ((*this)[i].grad) :
00610         (*this)[i].grad;
00611       for (size_t j = 0; j < group_for_fit->size(); j++) {
00612         group_for_fit->fit_gradients[j] +=
00613           (-1.0)/(cvm::real (group_for_fit->size())) *
00614           atom_grad;
00615       }
00616     }
00617   }
00618 
00619   if (b_rotate) {
00620     
00621     // add the rotation matrix contribution to the gradients
00622     cvm::rotation const rot_inv = rot.inverse();
00623     cvm::atom_pos const cog = this->center_of_geometry();
00624     
00625     for (size_t i = 0; i < this->size(); i++) {
00626 
00627       cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));
00628 
00629       for (size_t j = 0; j < group_for_fit->size(); j++) {
00630         // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
00631         cvm::quaternion const dxdq =
00632           rot.q.position_derivative_inner (pos_orig, (*this)[i].grad);
00633         // multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
00634         for (size_t iq = 0; iq < 4; iq++) {
00635           group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
00636         }
00637       }
00638     }
00639   }
00640   if (cvm::debug())
00641     cvm::log ("Done calculating fit gradients.\n");
00642 }
00643 
00644 
00645 std::vector<cvm::atom_pos> cvm::atom_group::positions() const
00646 {
00647   if (b_dummy)
00648     cvm::fatal_error ("Error: positions are not available "
00649                       "from a dummy atom group.\n");
00650 
00651   std::vector<cvm::atom_pos> x (this->size(), 0.0);
00652   cvm::atom_const_iter ai = this->begin();
00653   std::vector<cvm::atom_pos>::iterator xi = x.begin();
00654   for ( ; ai != this->end(); xi++, ai++) {
00655     *xi = ai->pos;
00656   }
00657   return x;
00658 }
00659 
00660 std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted (cvm::rvector const &shift) const
00661 {
00662   if (b_dummy)
00663     cvm::fatal_error ("Error: positions are not available "
00664                       "from a dummy atom group.\n");
00665 
00666   std::vector<cvm::atom_pos> x (this->size(), 0.0);
00667   cvm::atom_const_iter ai = this->begin();
00668   std::vector<cvm::atom_pos>::iterator xi = x.begin();
00669   for ( ; ai != this->end(); xi++, ai++) {
00670     *xi = (ai->pos + shift);
00671   }
00672   return x;
00673 }
00674 
00675 std::vector<cvm::rvector> cvm::atom_group::velocities() const
00676 {
00677   if (b_dummy)
00678     cvm::fatal_error ("Error: velocities are not available "
00679                       "from a dummy atom group.\n");
00680 
00681   std::vector<cvm::rvector> v (this->size(), 0.0);
00682   cvm::atom_const_iter ai = this->begin();
00683   std::vector<cvm::atom_pos>::iterator vi = v.begin();
00684   for ( ; ai != this->end(); vi++, ai++) {
00685     *vi = ai->vel;
00686   }
00687   return v;
00688 }
00689 
00690 std::vector<cvm::rvector> cvm::atom_group::system_forces() const
00691 {
00692   if (b_dummy)
00693     cvm::fatal_error ("Error: system forces are not available "
00694                       "from a dummy atom group.\n");
00695 
00696   std::vector<cvm::rvector> f (this->size(), 0.0);
00697   cvm::atom_const_iter ai = this->begin();
00698   std::vector<cvm::atom_pos>::iterator fi = f.begin();
00699   for ( ; ai != this->end(); fi++, ai++) {
00700     *fi = ai->system_force;
00701   }
00702   return f;
00703 }
00704 
00705 cvm::rvector cvm::atom_group::system_force() const
00706 {
00707   if (b_dummy)
00708     cvm::fatal_error ("Error: system forces are not available "
00709                       "from a dummy atom group.\n");
00710 
00711   cvm::rvector f (0.0);
00712   for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
00713     f += ai->system_force;
00714   }
00715   return f;
00716 }
00717 
00718 
00719 void cvm::atom_group::apply_colvar_force (cvm::real const &force)
00720 {
00721   if (b_dummy)
00722     return;
00723 
00724   if (noforce)
00725     cvm::fatal_error ("Error: sending a force to a group that has "
00726                       "\"enableForces\" set to off.\n");
00727 
00728   if (b_rotate) {
00729 
00730     // rotate forces back to the original frame
00731     cvm::rotation const rot_inv = rot.inverse();
00732     for (cvm::atom_iter ai = this->begin();
00733          ai != this->end(); ai++) {
00734       ai->apply_force (rot_inv.rotate (force * ai->grad));
00735     }
00736 
00737   } else {
00738 
00739     for (cvm::atom_iter ai = this->begin();
00740          ai != this->end(); ai++) {
00741       ai->apply_force (force * ai->grad);
00742     }
00743   }
00744 
00745   if ((b_center || b_rotate) && b_fit_gradients) {
00746 
00747     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
00748 
00749     // add the contribution from the roto-translational fit to the gradients
00750     if (b_rotate) {
00751       // rotate forces back to the original frame
00752       cvm::rotation const rot_inv = rot.inverse();
00753       for (size_t j = 0; j < group_for_fit->size(); j++) {
00754         (*group_for_fit)[j].apply_force (rot_inv.rotate (force * group_for_fit->fit_gradients[j]));
00755       }
00756     } else {
00757       for (size_t j = 0; j < group_for_fit->size(); j++) {
00758         (*group_for_fit)[j].apply_force (force * group_for_fit->fit_gradients[j]);
00759       }
00760     }
00761   }
00762 
00763 }
00764 
00765 
00766 void cvm::atom_group::apply_force (cvm::rvector const &force)
00767 {
00768   if (b_dummy)
00769     return;
00770 
00771   if (noforce)
00772     cvm::fatal_error ("Error: sending a force to a group that has "
00773                       "\"disableForces\" defined.\n");
00774 
00775   if (b_rotate) {
00776 
00777     cvm::rotation const rot_inv = rot.inverse();
00778     for (cvm::atom_iter ai = this->begin();
00779          ai != this->end(); ai++) {
00780       ai->apply_force (rot_inv.rotate ((ai->mass/this->total_mass) * force));
00781     }
00782 
00783   } else {
00784 
00785     for (cvm::atom_iter ai = this->begin();
00786          ai != this->end(); ai++) {
00787       ai->apply_force ((ai->mass/this->total_mass) * force);
00788     } 
00789   }
00790 }
00791 
00792 
00793 void cvm::atom_group::apply_forces (std::vector<cvm::rvector> const &forces)
00794 {
00795   if (b_dummy)
00796     return;
00797 
00798   if (noforce)
00799     cvm::fatal_error ("Error: sending a force to a group that has "
00800                       "\"disableForces\" defined.\n");
00801 
00802   if (forces.size() != this->size()) {
00803     cvm::fatal_error ("Error: trying to apply an array of forces to an atom "
00804                       "group which does not have the same length.\n");
00805   }
00806 
00807   if (b_rotate) {
00808 
00809     cvm::rotation const rot_inv = rot.inverse();
00810     cvm::atom_iter ai = this->begin();
00811     std::vector<cvm::rvector>::const_iterator fi = forces.begin();
00812     for ( ; ai != this->end(); fi++, ai++) {
00813       ai->apply_force (rot_inv.rotate (*fi));
00814     }
00815 
00816   } else {
00817 
00818     cvm::atom_iter ai = this->begin();
00819     std::vector<cvm::rvector>::const_iterator fi = forces.begin();
00820     for ( ; ai != this->end(); fi++, ai++) {
00821       ai->apply_force (*fi);
00822     }
00823   }
00824 }
00825 

Generated on Sat May 25 04:07:14 2013 for NAMD by  doxygen 1.3.9.1