Difference for src/colvaratoms.C from version 1.19 to 1.20

version 1.19version 1.20
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include "colvarmodule.h" #include "colvarmodule.h"
 #include "colvarparse.h" #include "colvarparse.h"
 #include "colvaratoms.h" #include "colvaratoms.h"
  
  cvm::atom::atom()
  {
    index = -1;
    id = -1;
    reset_data();
  }
  
  
  cvm::atom::atom(int atom_number)
  {
    colvarproxy *p = cvm::proxy;
    index = p->init_atom(atom_number);
    if (cvm::debug()) {
      cvm::log("The index of this atom in the colvarproxy arrays is "+
               cvm::to_str(index)+".\n");
    }
    id = p->get_atom_id(index);
    update_mass();
    reset_data();
  }
  
  
  cvm::atom::atom(cvm::residue_id const &residue,
                  std::string const     &atom_name,
                  std::string const     &segment_id)
  {
    colvarproxy *p = cvm::proxy;
    index = p->init_atom(residue, atom_name, segment_id);
    if (cvm::debug()) {
      cvm::log("The index of this atom in the colvarproxy_namd arrays is "+
               cvm::to_str(index)+".\n");
    }
    id = p->get_atom_id(index);
    update_mass();
    reset_data();
  }
  
  
  cvm::atom::atom(atom const &a)
    : index(a.index)
  {
    id = (cvm::proxy)->get_atom_id(index);
    update_mass();
    reset_data();
  }
  
  
 // member functions of the "atom" class depend tightly on the MD interface, and are cvm::atom::~atom()
 // thus defined in colvarproxy_xxx.cpp {
    if (index >= 0) {
      (cvm::proxy)->clear_atom(index);
    }
  }
  
 // in this file only atom_group functions are defined 
  
  
  // TODO change this arrangement
 // Note: "conf" is the configuration of the cvc who is using this atom group; // Note: "conf" is the configuration of the cvc who is using this atom group;
 // "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...) // "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
 cvm::atom_group::atom_group(std::string const &conf, cvm::atom_group::atom_group(std::string const &conf,
                              char const        *key)                             char const        *key_in)
   : b_center(false), b_rotate(false), b_user_defined_fit(false), 
     b_fit_gradients(false), 
     ref_pos_group(NULL), 
     noforce(false) 
 { {
   cvm::log("Defining atom group \""+   key = key_in;
             std::string(key)+"\".\n");   cvm::log("Defining atom group \"" + key + "\".\n");
    init();
   // real work is done by parse   // real work is done by parse
   parse(conf, key);   parse(conf);
    setup();
 } }
  
  
 cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms) cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms_in)
   : b_dummy(false), b_center(false), b_rotate(false), 
     b_fit_gradients(false), ref_pos_group(NULL), 
     noforce(false) 
 { {
   this->reserve(atoms.size());   init();
   for (size_t i = 0; i < atoms.size(); i++) {   atoms = atoms_in;
     this->push_back(atoms[i]);   setup();
   } 
   total_mass = 0.0; 
   for (cvm::atom_iter ai = this->begin(); 
        ai != this->end(); ai++) { 
     total_mass += ai->mass; 
   } 
 } }
  
  
 cvm::atom_group::atom_group() cvm::atom_group::atom_group()
   : b_dummy(false), b_center(false), b_rotate(false), 
     b_user_defined_fit(false), b_fit_gradients(false), 
     ref_pos_group(NULL), noforce(false) 
 { {
   total_mass = 0.0;   init();
 } }
  
  
 cvm::atom_group::~atom_group() cvm::atom_group::~atom_group()
 { {
   if (ref_pos_group) {   if (is_enabled(f_ag_scalable)) {
     delete ref_pos_group;     (cvm::proxy)->clear_atom_group(index);
     ref_pos_group = NULL;     index = -1;
    }
  
    if (fitting_group) {
      delete fitting_group;
      fitting_group = NULL;
   }   }
 } }
  
  
 void cvm::atom_group::add_atom(cvm::atom const &a) int cvm::atom_group::add_atom(cvm::atom const &a)
 { {
   if (b_dummy) {   if (a.id < 0) {
     cvm::error("Error: cannot add atoms to a dummy group.\n", INPUT_ERROR);     return COLVARS_ERROR;
   } else {   }
     this->push_back(a); 
    for (size_t i = 0; i < atoms_ids.size(); i++) {
      if (atoms_ids[i] == a.id) {
        if (cvm::debug())
          cvm::log("Discarding doubly counted atom with number "+
                   cvm::to_str(a.id+1)+".\n");
        return COLVARS_OK;
      }
    }
  
    // for consistency with add_atom_id(), we update the list as well
    atoms_ids.push_back(a.id);
    atoms.push_back(a);
     total_mass += a.mass;     total_mass += a.mass;
    total_charge += a.charge;
  
    return COLVARS_OK;
  }
  
  
  int cvm::atom_group::add_atom_id(int aid)
  {
    if (aid < 0) {
      return COLVARS_ERROR;
    }
  
    for (size_t i = 0; i < atoms_ids.size(); i++) {
      if (atoms_ids[i] == aid) {
        if (cvm::debug())
          cvm::log("Discarding doubly counted atom with number "+
                   cvm::to_str(aid+1)+".\n");
        return COLVARS_OK;
   }   }
 } }
  
    atoms_ids.push_back(aid);
    return COLVARS_OK;
  }
  
 void cvm::atom_group::reset_mass(std::string &name, int i, int j) 
  int cvm::atom_group::remove_atom(cvm::atom_iter ai)
  {
    if (is_enabled(f_ag_scalable)) {
      cvm::error("Error: cannot remove atoms from a scalable group.\n", INPUT_ERROR);
      return COLVARS_ERROR;
    }
  
    if (!this->size()) {
      cvm::error("Error: trying to remove an atom from an empty group.\n", INPUT_ERROR);
      return COLVARS_ERROR;
    } else {
      total_mass -= ai->mass;
      total_charge -= ai->charge;
      atoms_ids.erase(atoms_ids.begin() + (ai - atoms.begin()));
      atoms.erase(ai);
    }
  
    return COLVARS_OK;
  }
  
  
  int cvm::atom_group::init()
  {
    if (!key.size()) key = "atoms";
    atoms.clear();
  
    // TODO: check with proxy whether atom forces etc are available
    init_ag_requires();
  
    index = -1;
  
    b_center = false;
    b_rotate = false;
    b_user_defined_fit = false;
    b_fit_gradients = false;
    fitting_group = NULL;
  
    noforce = false;
  
    total_mass = 0.0;
    total_charge = 0.0;
  
    cog.reset();
    com.reset();
  
    return COLVARS_OK;
  }
  
  
  int cvm::atom_group::setup()
 { {
    for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
      ai->update_mass();
      ai->update_charge();
    }
    update_total_mass();
    update_total_charge();
    return COLVARS_OK;
  }
  
  
  void cvm::atom_group::update_total_mass()
  {
    if (b_dummy) {
      total_mass = 1.0;
      return;
    }
  
    if (is_enabled(f_ag_scalable)) {
      total_mass = (cvm::proxy)->get_atom_group_mass(index);
    } else {
   total_mass = 0.0;   total_mass = 0.0;
   for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
        ai != this->end(); ai++) { 
     total_mass += ai->mass;     total_mass += ai->mass;
   }   }
    }
  }
  
  
  void cvm::atom_group::reset_mass(std::string &name, int i, int j)
  {
    update_total_mass();
   cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+   cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+
             cvm::to_str(j)+". "+ cvm::to_str(this->size())+            cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+
             " atoms: total mass = "+cvm::to_str(this->total_mass)+".\n");            " atoms: total mass = "+cvm::to_str(total_mass)+".\n");
  }
  
  
  void cvm::atom_group::update_total_charge()
  {
    if (b_dummy) {
      total_charge = 0.0;
      return;
    }
  
    if (is_enabled(f_ag_scalable)) {
      total_charge = (cvm::proxy)->get_atom_group_charge(index);
    } else {
      total_charge = 0.0;
      for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
        total_charge += ai->charge;
      }
    }
 } }
  
 int cvm::atom_group::parse(std::string const &conf, 
                              char const        *key) int cvm::atom_group::parse(std::string const &conf)
 { {
   std::string group_conf;   std::string group_conf;
  
    // TODO move this to the cvc class constructor/init
  
   // save_delimiters is set to false for this call, because "conf" is   // save_delimiters is set to false for this call, because "conf" is
   // not the config string of this group, but of its parent object   // not the config string of this group, but of its parent object
   // (which has already taken care of the delimiters)   // (which has already taken care of the delimiters)
   save_delimiters = false;   save_delimiters = false;
   key_lookup(conf, key, group_conf, dummy_pos);   key_lookup(conf, key.c_str(), group_conf, dummy_pos);
   // restoring the normal value, because we do want keywords checked   // restoring the normal value, because we do want keywords checked
   // inside "group_conf"   // inside "group_conf"
   save_delimiters = true;   save_delimiters = true;
  
   if (group_conf.size() == 0) {   if (group_conf.size() == 0) {
     cvm::error("Error: atom group \""+std::string(key)+     cvm::error("Error: atom group \""+key+
                 "\" is set, but has no definition.\n",                 "\" is set, but has no definition.\n",
                 INPUT_ERROR);                 INPUT_ERROR);
     return COLVARS_ERROR;     return COLVARS_ERROR;
Line 108
Line 278
  
   cvm::increase_depth();   cvm::increase_depth();
  
   cvm::log("Initializing atom group \""+std::string(key)+"\".\n");   cvm::log("Initializing atom group \""+key+"\".\n");
  
    description = "atom group " + key;
  
   // whether or not to include messages in the log   // whether or not to include messages in the log
   // colvarparse::Parse_Mode mode = parse_silent;   // colvarparse::Parse_Mode mode = parse_silent;
Line 117
Line 289
   //   get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);   //   get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
   //   if (b_verbose) mode = parse_normal;   //   if (b_verbose) mode = parse_normal;
   // }   // }
   colvarparse::Parse_Mode mode = parse_normal;   // colvarparse::Parse_Mode mode = parse_normal;
  
    int parse_error = COLVARS_OK;
  
   {   {
     //    std::vector<int> atom_indexes; 
     std::string numbers_conf = "";     std::string numbers_conf = "";
     size_t pos = 0;     size_t pos = 0;
     std::vector<int> atom_indexes; 
     while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {     while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
        parse_error |= add_atom_numbers(numbers_conf);
        numbers_conf = "";
      }
    }
  
    {
      std::string index_group_name;
      if (get_keyval(group_conf, "indexGroup", index_group_name)) {
        // use an index group from the index file read globally
        parse_error |= add_index_group(index_group_name);
      }
    }
  
    {
      std::string range_conf = "";
      size_t pos = 0;
      while (key_lookup(group_conf, "atomNumbersRange",
                        range_conf, pos)) {
        parse_error |= add_atom_numbers_range(range_conf);
        range_conf = "";
      }
    }
  
    {
      std::vector<std::string> psf_segids;
      get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string>());
      std::vector<std::string>::iterator psii;
      for (psii = psf_segids.begin(); psii < psf_segids.end(); ++psii) {
        if ( (psii->size() == 0) || (psii->size() > 4) ) {
          cvm::error("Error: invalid PSF segment identifier provided, \""+
                     (*psii)+"\".\n", INPUT_ERROR);
        }
      }
  
      std::string range_conf = "";
      size_t pos = 0;
      size_t range_count = 0;
      psii = psf_segids.begin();
      while (key_lookup(group_conf, "atomNameResidueRange",
                        range_conf, pos)) {
        range_count++;
        if (psf_segids.size() && (range_count > psf_segids.size())) {
          cvm::error("Error: more instances of \"atomNameResidueRange\" than "
                     "values of \"psfSegID\".\n", INPUT_ERROR);
        } else {
          parse_error |= add_atom_name_residue_range(psf_segids.size() ?
            *psii : std::string(""), range_conf);
          if (psf_segids.size()) psii++;
        }
        range_conf = "";
      }
    }
  
    {
      // read the atoms from a file
      std::string atoms_file_name;
      if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""))) {
  
        std::string atoms_col;
        if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""))) {
          cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n",
                     INPUT_ERROR);
        }
  
        double atoms_col_value;
        bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0);
        if (atoms_col_value_defined && (!atoms_col_value)) {
          cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
        }
  
        cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
      }
    }
  
    // Catch any errors from all the initialization steps above
    if (parse_error || cvm::get_error()) return (parse_error || cvm::get_error());
  
    // checks of doubly-counted atoms have been handled by add_atom() already
  
    if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos())) {
      b_dummy = true;
      // note: atoms_ids.size() is used here in lieu of atoms.size(),
      // which can be empty for scalable groups
      if (atoms_ids.size()) {
        cvm::error("Error: cannot set up group \""+
                   key+"\" as a dummy atom "
                   "and provide it with atom definitions.\n", INPUT_ERROR);
      }
    } else {
      b_dummy = false;
  
      if (!(atoms_ids.size())) {
        cvm::error("Error: no atoms defined for atom group \""+
                   key+"\".\n", INPUT_ERROR);
      }
  
      // whether these atoms will ever receive forces or not
      bool enable_forces = true;
      // disableForces is deprecated
      if (get_keyval(group_conf, "enableForces", enable_forces, true)) {
        noforce = !enable_forces;
      } else {
        get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
      }
    }
  
    if (is_enabled(f_ag_scalable) && !b_dummy) {
      index = (cvm::proxy)->init_atom_group(atoms_ids);
    }
  
    parse_error |= parse_fitting_options(group_conf);
  
    // TODO move this to colvarparse object
    check_keywords(group_conf, key.c_str());
    if (cvm::get_error()) {
      cvm::error("Error setting up atom group \""+key+"\".");
      return COLVARS_ERROR;
    }
  
    // Calculate all required properties (such as total mass)
    setup();
  
    if (cvm::debug())
      cvm::log("Done initializing atom group \""+key+"\".\n");
  
    cvm::log("Atom group \""+key+"\" defined, "+
              cvm::to_str(atoms_ids.size())+" atoms initialized: total mass = "+
       cvm::to_str(total_mass)+", total charge = "+
              cvm::to_str(total_charge)+".\n");
  
    cvm::decrease_depth();
  
    return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
  }
  
  
  int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
  {
    std::vector<int> atom_indexes;
  
       if (numbers_conf.size()) {       if (numbers_conf.size()) {
         std::istringstream is(numbers_conf);         std::istringstream is(numbers_conf);
         int ia;         int ia;
Line 134
Line 446
       }       }
  
       if (atom_indexes.size()) {       if (atom_indexes.size()) {
         this->reserve(this->size()+atom_indexes.size());     atoms_ids.reserve(atoms_ids.size()+atom_indexes.size());
  
      if (is_enabled(f_ag_scalable)) {
        for (size_t i = 0; i < atom_indexes.size(); i++) {
          add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i]));
        }
      } else {
        // if we are handling the group on rank 0, better allocate the vector in one shot
        atoms.reserve(atoms.size()+atom_indexes.size());
         for (size_t i = 0; i < atom_indexes.size(); i++) {         for (size_t i = 0; i < atom_indexes.size(); i++) {
           this->push_back(cvm::atom(atom_indexes[i]));         add_atom(cvm::atom(atom_indexes[i]));
         }         }
      }
  
         if (cvm::get_error()) return COLVARS_ERROR;         if (cvm::get_error()) return COLVARS_ERROR;
       } else {       } else {
         cvm::error("Error: no numbers provided for \""         cvm::error("Error: no numbers provided for \""
Line 145
Line 467
         return COLVARS_ERROR;         return COLVARS_ERROR;
       }       }
  
       atom_indexes.clear();   return COLVARS_OK;
     }     }
  
     std::string index_group_name; 
     if (get_keyval(group_conf, "indexGroup", index_group_name)) { int cvm::atom_group::add_index_group(std::string const &index_group_name)
       // use an index group from the index file read globally {
       std::list<std::string>::iterator names_i = cvm::index_group_names.begin();       std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
       std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();       std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
       for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) {       for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) {
         if (*names_i == index_group_name)         if (*names_i == index_group_name)
           break;           break;
       }       }
  
       if (names_i == cvm::index_group_names.end()) {       if (names_i == cvm::index_group_names.end()) {
         cvm::error("Error: could not find index group "+         cvm::error("Error: could not find index group "+
                     index_group_name+" among those provided by the index file.\n",                     index_group_name+" among those provided by the index file.\n",
                     INPUT_ERROR);                     INPUT_ERROR);
         return COLVARS_ERROR;         return COLVARS_ERROR;
       }       }
       this->reserve(index_groups_i->size()); 
    atoms_ids.reserve(atoms_ids.size()+index_groups_i->size());
  
    if (is_enabled(f_ag_scalable)) {
       for (size_t i = 0; i < index_groups_i->size(); i++) {       for (size_t i = 0; i < index_groups_i->size(); i++) {
         this->push_back(cvm::atom((*index_groups_i)[i]));       add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i]));
       }       }
       if (cvm::get_error()) return COLVARS_ERROR;   } else {
      atoms.reserve(atoms.size()+index_groups_i->size());
      for (size_t i = 0; i < index_groups_i->size(); i++) {
        add_atom(cvm::atom((*index_groups_i)[i]));
     }     }
   }   }
  
   {   if (cvm::get_error())
     std::string range_conf = "";     return COLVARS_ERROR;
     size_t pos = 0; 
     while (key_lookup(group_conf, "atomNumbersRange",   return COLVARS_OK;
                        range_conf, pos)) { }
  
  
  int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf)
  {
       if (range_conf.size()) {       if (range_conf.size()) {
         std::istringstream is(range_conf);         std::istringstream is(range_conf);
         int initial, final;         int initial, final;
Line 184
Line 516
         if ( (is >> initial) && (initial > 0) &&         if ( (is >> initial) && (initial > 0) &&
              (is >> dash) && (dash == '-') &&              (is >> dash) && (dash == '-') &&
              (is >> final) && (final > 0) ) {              (is >> final) && (final > 0) ) {
  
        atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
  
        if (is_enabled(f_ag_scalable)) {
           for (int anum = initial; anum <= final; anum++) {           for (int anum = initial; anum <= final; anum++) {
             this->push_back(cvm::atom(anum));           add_atom_id((cvm::proxy)->check_atom_id(anum));
           }           }
           if (cvm::get_error()) return COLVARS_ERROR;       } else {
           range_conf = "";         atoms.reserve(atoms.size() + (final - initial + 1));
           continue;         for (int anum = initial; anum <= final; anum++) {
            add_atom(cvm::atom(anum));
         }         }
       }       }
  
      }
      if (cvm::get_error()) return COLVARS_ERROR;
    } else {
       cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+       cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+
                     range_conf+"\".\n", INPUT_ERROR);                     range_conf+"\".\n", INPUT_ERROR);
     }     return COLVARS_ERROR;
   } 
  
   std::vector<std::string> psf_segids; 
   get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode); 
   for (std::vector<std::string>::iterator psii = psf_segids.begin(); 
        psii < psf_segids.end(); ++psii) { 
  
     if ( (psii->size() == 0) || (psii->size() > 4) ) { 
       cvm::error("Error: invalid segmend identifier provided, \""+ 
                   (*psii)+"\".\n", INPUT_ERROR); 
     } 
   }   }
  
   {   return COLVARS_OK;
     std::string range_conf = ""; 
     size_t pos = 0; 
     size_t range_count = 0; 
     std::vector<std::string>::iterator psii = psf_segids.begin(); 
     while (key_lookup(group_conf, "atomNameResidueRange", 
                        range_conf, pos)) { 
       range_count++; 
  
       if (range_count > psf_segids.size()) { 
         cvm::error("Error: more instances of \"atomNameResidueRange\" than " 
                     "values of \"psfSegID\".\n", INPUT_ERROR); 
       }       }
  
       std::string const &psf_segid = psf_segids.size() ? *psii : std::string(""); 
  
  int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
                                                   std::string const &range_conf)
  {
       if (range_conf.size()) {       if (range_conf.size()) {
  
         std::istringstream is(range_conf);         std::istringstream is(range_conf);
         std::string atom_name;         std::string atom_name;
         int initial, final;         int initial, final;
Line 235
Line 554
              (is >> initial) && (initial > 0) &&              (is >> initial) && (initial > 0) &&
              (is >> dash) && (dash == '-') &&              (is >> dash) && (dash == '-') &&
              (is >> final) && (final > 0) ) {              (is >> final) && (final > 0) ) {
  
        atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
  
        if (is_enabled(f_ag_scalable)) {
           for (int resid = initial; resid <= final; resid++) {           for (int resid = initial; resid <= final; resid++) {
             this->push_back(cvm::atom(resid, atom_name, psf_segid));           add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid));
          }
        } else {
          atoms.reserve(atoms.size() + (final - initial + 1));
          for (int resid = initial; resid <= final; resid++) {
            add_atom(cvm::atom(resid, atom_name, psf_segid));
          }
           }           }
  
           if (cvm::get_error()) return COLVARS_ERROR;           if (cvm::get_error()) return COLVARS_ERROR;
           range_conf = ""; 
         } else {         } else {
           cvm::error("Error: cannot parse definition for \""           cvm::error("Error: cannot parse definition for \""
                       "atomNameResidueRange\", \""+                       "atomNameResidueRange\", \""+
                       range_conf+"\".\n");                       range_conf+"\".\n");
        return COLVARS_ERROR;
         }         }
  
       } else {       } else {
         cvm::error("Error: atomNameResidueRange with empty definition.\n");         cvm::error("Error: atomNameResidueRange with empty definition.\n");
      return COLVARS_ERROR;
       }       }
    return COLVARS_OK;
       if (psf_segid.size()) 
         ++psii; 
     } 
   } 
  
   { 
     // read the atoms from a file 
     std::string atoms_file_name; 
     if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""), mode)) { 
  
       std::string atoms_col; 
       if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""), mode)) { 
         cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n", 
                     INPUT_ERROR); 
       } 
  
       double atoms_col_value; 
       bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0, mode); 
       if (atoms_col_value_defined && (!atoms_col_value)) { 
         cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR); 
       } 
  
       cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value); 
     } 
   } 
  
   // Catch any errors from all the initialization steps above 
   if (cvm::get_error()) return COLVARS_ERROR; 
  
   for (std::vector<cvm::atom>::iterator a1 = this->begin(); 
        a1 != this->end(); ++a1) { 
     std::vector<cvm::atom>::iterator a2 = a1; 
     ++a2; 
     for ( ; a2 != this->end(); ++a2) { 
       if (a1->id == a2->id) { 
         if (cvm::debug()) 
           cvm::log("Discarding doubly counted atom with number "+ 
                     cvm::to_str(a1->id+1)+".\n"); 
         a2 = this->erase(a2); 
         if (a2 == this->end()) 
           break; 
       } 
     } 
   } 
  
   if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) { 
     b_dummy = true; 
     this->total_mass = 1.0; 
   } else 
     b_dummy = false; 
  
   if (b_dummy && (this->size())) { 
     cvm::error("Error: cannot set up group \""+ 
                 std::string(key)+"\" as a dummy atom " 
                 "and provide it with atom definitions.\n", INPUT_ERROR); 
   } 
  
 #if(! defined(COLVARS_STANDALONE)) 
   if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) { 
     cvm::error("Error: no atoms defined for atom group \""+ 
                       std::string(key)+"\".\n"); 
   } 
 #endif 
  
   if (!b_dummy) { 
  
     // calculate total mass (TODO: this is the step that most needs deferred re-initialization) 
     this->total_mass = 0.0; 
     for (cvm::atom_iter ai = this->begin(); 
          ai != this->end(); ai++) { 
       this->total_mass += ai->mass; 
     } 
  
     // whether these atoms will ever receive forces or not 
     bool enable_forces = true; 
     // disableForces is deprecated 
     if (get_keyval(group_conf, "enableForces", enable_forces, true, mode)) { 
       noforce = !enable_forces; 
     } else { 
       get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent); 
     } 
   }   }
  
   // FITTING OPTIONS 
  
   bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false, mode); int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
   bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false, mode); {
    bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
    bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false);
   // is the user setting explicit options?   // is the user setting explicit options?
   b_user_defined_fit = b_defined_center || b_defined_rotate;   b_user_defined_fit = b_defined_center || b_defined_rotate;
  
   get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true, mode);   get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true);
  
   if (b_center || b_rotate) {   if (b_center || b_rotate) {
  
Line 348
Line 598
       cvm::error("Error: centerReference or rotateReference "       cvm::error("Error: centerReference or rotateReference "
                   "cannot be defined for a dummy atom.\n");                   "cannot be defined for a dummy atom.\n");
  
      bool b_ref_pos_group = false;
     if (key_lookup(group_conf, "refPositionsGroup")) {     if (key_lookup(group_conf, "refPositionsGroup")) {
        b_ref_pos_group = true;
        cvm::log("Warning: keyword \"refPositionsGroup\" is deprecated: please use \"fittingGroup\" instead.\n");
      }
  
      if (b_ref_pos_group || key_lookup(group_conf, "fittingGroup")) {
       // instead of this group, define another group to compute the fit       // instead of this group, define another group to compute the fit
       if (ref_pos_group) {       if (fitting_group) {
         cvm::error("Error: the atom group \""+         cvm::error("Error: the atom group \""+
                     std::string(key)+"\" has already a reference group "                    key+"\" has already a reference group "
                     "for the rototranslational fit, which was communicated by the "                     "for the rototranslational fit, which was communicated by the "
                     "colvar component.  You should not use refPositionsGroup "                    "colvar component.  You should not use fittingGroup "
                     "in this case.\n");                     "in this case.\n");
       }       }
       cvm::log("Within atom group \""+std::string(key)+"\":\n");       cvm::log("Within atom group \""+key+"\":\n");
       ref_pos_group = new atom_group(group_conf, "refPositionsGroup");       fitting_group = b_ref_pos_group ?
          new atom_group(group_conf, "refPositionsGroup") :
          new atom_group(group_conf, "fittingGroup");
  
       // regardless of the configuration, fit gradients must be calculated by refPositionsGroup       // regardless of the configuration, fit gradients must be calculated by fittingGroup
       ref_pos_group->b_fit_gradients = this->b_fit_gradients;       fitting_group->b_fit_gradients = this->b_fit_gradients;
       this->b_fit_gradients = false; 
     }     }
  
     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;     atom_group *group_for_fit = fitting_group ? fitting_group : this;
  
     get_keyval(group_conf, "refPositions", ref_pos, ref_pos, mode);     get_keyval(group_conf, "refPositions", ref_pos, ref_pos);
  
     std::string ref_pos_file;     std::string ref_pos_file;
     if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""), mode)) {     if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""))) {
  
       if (ref_pos.size()) {       if (ref_pos.size()) {
         cvm::error("Error: cannot specify \"refPositionsFile\" and "         cvm::error("Error: cannot specify \"refPositionsFile\" and "
Line 380
Line 637
       std::string ref_pos_col;       std::string ref_pos_col;
       double ref_pos_col_value=0.0;       double ref_pos_col_value=0.0;
  
       if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""), mode)) {       if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""))) {
         // if provided, use PDB column to select coordinates         // if provided, use PDB column to select coordinates
         bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);         bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0);
         if (found && ref_pos_col_value == 0.0)         if (found && ref_pos_col_value == 0.0) {
           cvm::error("Error: refPositionsColValue, "           cvm::error("Error: refPositionsColValue, "
                       "if provided, must be non-zero.\n");                      "if provided, must be non-zero.\n", INPUT_ERROR);
            return COLVARS_ERROR;
          }
       } else {       } else {
         // if not, rely on existing atom indices for the group         // if not, rely on existing atom indices for the group
         group_for_fit->create_sorted_ids();         group_for_fit->create_sorted_ids();
Line 403
Line 662
           cvm::error("Error: the number of reference positions provided("+           cvm::error("Error: the number of reference positions provided("+
                       cvm::to_str(ref_pos.size())+                       cvm::to_str(ref_pos.size())+
                       ") does not match the number of atoms within \""+                       ") does not match the number of atoms within \""+
                       std::string(key)+                      key+
                       "\" ("+cvm::to_str(group_for_fit->size())+                       "\" ("+cvm::to_str(group_for_fit->size())+
                       "): to perform a rotational fit, "+                       "): to perform a rotational fit, "+
                       "these numbers should be equal.\n", INPUT_ERROR);                       "these numbers should be equal.\n", INPUT_ERROR);
Line 413
Line 672
       center_ref_pos();       center_ref_pos();
  
     } else {     } else {
 #if(! defined(COLVARS_STANDALONE))       cvm::error("Error: no reference positions provided.\n", INPUT_ERROR);
       cvm::error("Error: no reference positions provided.\n");       return COLVARS_ERROR;
 #endif 
     }     }
  
     if (b_fit_gradients) {     if (b_fit_gradients) {
Line 424
Line 682
     }     }
  
     if (b_rotate && !noforce) {     if (b_rotate && !noforce) {
       cvm::log("Warning: atom group \""+std::string(key)+       cvm::log("Warning: atom group \""+key+
                 "\" will be aligned to a fixed orientation given by the reference positions provided.  "                 "\" will be aligned to a fixed orientation given by the reference positions provided.  "
                 "If the internal structure of the group changes too much (i.e. its RMSD is comparable "                 "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
                 "to its radius of gyration), the optimal rotation and its gradients may become discontinuous.  "                 "to its radius of gyration), the optimal rotation and its gradients may become discontinuous.  "
                 "If that happens, use refPositionsGroup (or a different definition for it if already defined) "                "If that happens, use fittingGroup (or a different definition for it if already defined) "
                 "to align the coordinates.\n");                 "to align the coordinates.\n");
       // initialize rot member data       // initialize rot member data
       rot.request_group1_gradients(this->size());       rot.request_group1_gradients(group_for_fit->size());
     } 
  
   }   }
  
   if (cvm::debug()) 
     cvm::log("Done initializing atom group with name \""+ 
               std::string(key)+"\".\n"); 
  
   this->check_keywords(group_conf, key); 
   if (cvm::get_error()) { 
     cvm::error("Error setting up atom group \""+std::string(key)+"\"."); 
     return COLVARS_ERROR; 
   }   }
  
   cvm::log("Atom group \""+std::string(key)+"\" defined, "+   return COLVARS_OK;
             cvm::to_str(this->size())+" atoms initialized: total mass = "+ 
             cvm::to_str(this->total_mass)+".\n"); 
  
   cvm::decrease_depth(); 
  
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); 
 } }
  
  
Line 463
Line 704
     return COLVARS_OK;     return COLVARS_OK;
  
   std::list<int> temp_id_list;   std::list<int> temp_id_list;
   cvm::atom_iter ai;   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
   for (ai = this->begin(); ai != this->end(); ai++) { 
     temp_id_list.push_back(ai->id);     temp_id_list.push_back(ai->id);
   }   }
   temp_id_list.sort();   temp_id_list.sort();
Line 486
Line 726
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 } }
  
  
 void cvm::atom_group::center_ref_pos() void cvm::atom_group::center_ref_pos()
 { {
   // save the center of geometry of ref_pos and then subtract it from 
   // them; in this way it will be possible to use ref_pos also for 
   // the rotational fit 
   // This is called either by atom_group::parse or by CVCs that set 
   // reference positions (eg. RMSD, eigenvector) 
   ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);   ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);
   std::vector<cvm::atom_pos>::iterator pi;   std::vector<cvm::atom_pos>::iterator pi;
   for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {   for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
Line 504
Line 740
   }   }
 } }
  
  
 void cvm::atom_group::read_positions() void cvm::atom_group::read_positions()
 { {
   if (b_dummy) return;   if (b_dummy) return;
  
   for (cvm::atom_iter ai = this->begin();   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
        ai != this->end(); ai++) { 
     ai->read_position();     ai->read_position();
   }   }
  
   if (ref_pos_group)   if (fitting_group)
     ref_pos_group->read_positions();     fitting_group->read_positions();
 } }
  
  
  int cvm::atom_group::calc_required_properties()
  {
    // TODO check if the com is needed?
    calc_center_of_mass();
    calc_center_of_geometry();
  
    if (!is_enabled(f_ag_scalable)) {
      if (b_center || b_rotate) {
        if (fitting_group) {
          fitting_group->calc_center_of_geometry();
        }
  
        calc_apply_roto_translation();
  
        // update COM and COG after fitting
        calc_center_of_geometry();
        calc_center_of_mass();
        if (fitting_group) {
          fitting_group->calc_center_of_geometry();
        }
      }
    }
  
    // TODO calculate elements of scalable cvc's here before reduction
  
    return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
  }
  
  
 void cvm::atom_group::calc_apply_roto_translation() void cvm::atom_group::calc_apply_roto_translation()
 { {
   atom_group *fit_group = ref_pos_group ? ref_pos_group : this;   // store the laborarory-frame COGs for when they are needed later
    cog_orig = this->center_of_geometry();
    if (fitting_group) {
      fitting_group->cog_orig = fitting_group->center_of_geometry();
    }
  
   if (b_center) {   if (b_center) {
     // center on the origin first     // center on the origin first
     cvm::atom_pos const cog = fit_group->center_of_geometry();     cvm::atom_pos const rpg_cog = fitting_group ?
     for (cvm::atom_iter ai = this->begin();       fitting_group->center_of_geometry() : this->center_of_geometry();
          ai != this->end(); ai++) {     apply_translation(-1.0 * rpg_cog);
       ai->pos -= cog;     if (fitting_group) {
        fitting_group->apply_translation(-1.0 * rpg_cog);
     }     }
   }   }
  
   if (b_rotate) {   if (b_rotate) {
     // rotate the group (around the center of geometry if b_center is     // rotate the group (around the center of geometry if b_center is
     // true, around the origin otherwise)     // true, around the origin otherwise)
     rot.calc_optimal_rotation(fit_group->positions(), ref_pos);     rot.calc_optimal_rotation(fitting_group ?
                                fitting_group->positions() :
                                this->positions(),
                                ref_pos);
  
     for (cvm::atom_iter ai = this->begin();     cvm::atom_iter ai;
          ai != this->end(); ai++) {     for (ai = this->begin(); ai != this->end(); ai++) {
        ai->pos = rot.rotate(ai->pos);
      }
      if (fitting_group) {
        for (ai = fitting_group->begin(); ai != fitting_group->end(); ai++) {
       ai->pos = rot.rotate(ai->pos);       ai->pos = rot.rotate(ai->pos);
     }     }
   }   }
    }
  
   if (b_center) {   if (b_center) {
     // align with the center of geometry of ref_pos     // align with the center of geometry of ref_pos
     for (cvm::atom_iter ai = this->begin();     apply_translation(ref_pos_cog);
          ai != this->end(); ai++) {     if (fitting_group) {
       ai->pos += ref_pos_cog;       fitting_group->apply_translation(ref_pos_cog);
     }     }
   }   }
    // update of COM and COG is done from the calling routine
 } }
  
  
 void cvm::atom_group::apply_translation(cvm::rvector const &t) void cvm::atom_group::apply_translation(cvm::rvector const &t)
 { {
   if (b_dummy) return;   if (b_dummy) {
      cvm::error("Error: cannot translate the coordinates of a dummy atom group.\n", INPUT_ERROR);
   for (cvm::atom_iter ai = this->begin();     return;
        ai != this->end(); ai++) { 
     ai->pos += t; 
   } 
 } }
  
 void cvm::atom_group::apply_rotation(cvm::rotation const &rot)   if (is_enabled(f_ag_scalable)) {
 {     cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", INPUT_ERROR);
   if (b_dummy) return;     return;
    }
  
   for (cvm::atom_iter ai = this->begin();   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
        ai != this->end(); ai++) {     ai->pos += t;
     ai->pos = rot.rotate(ai->pos); 
   }   }
 } }
  
Line 577
Line 855
  
   if (b_rotate) {   if (b_rotate) {
  
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) { 
       ai->read_velocity();       ai->read_velocity();
       ai->vel = rot.rotate(ai->vel);       ai->vel = rot.rotate(ai->vel);
     }     }
  
   } else {   } else {
  
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) { 
       ai->read_velocity();       ai->read_velocity();
     }     }
   }   }
 } }
  
 void cvm::atom_group::read_system_forces() 
  // TODO make this a calc function
  void cvm::atom_group::read_total_forces()
 { {
   if (b_dummy) return;   if (b_dummy) return;
  
   if (b_rotate) {   if (b_rotate) {
  
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) {       ai->read_total_force();
       ai->read_system_force();       ai->total_force = rot.rotate(ai->total_force);
       ai->system_force = rot.rotate(ai->system_force); 
     }     }
  
   } else {   } else {
  
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) {       ai->read_total_force();
       ai->read_system_force(); 
     }     }
   }   }
 } }
  
 cvm::atom_pos cvm::atom_group::center_of_geometry() const 
 { 
   if (b_dummy) 
     return dummy_atom_pos; 
  
   cvm::atom_pos cog(0.0, 0.0, 0.0); int cvm::atom_group::calc_center_of_geometry()
   for (cvm::atom_const_iter ai = this->begin(); {
        ai != this->end(); ai++) {   if (b_dummy) {
      cog = dummy_atom_pos;
    } else {
      cog.reset();
      for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
     cog += ai->pos;     cog += ai->pos;
   }   }
   cog /= this->size();   cog /= this->size();
   return cog;   }
    return COLVARS_OK;
 } }
  
 cvm::atom_pos cvm::atom_group::center_of_mass() const 
 { 
   if (b_dummy) 
     return dummy_atom_pos; 
  
   cvm::atom_pos com(0.0, 0.0, 0.0); int cvm::atom_group::calc_center_of_mass()
   for (cvm::atom_const_iter ai = this->begin(); {
        ai != this->end(); ai++) {   if (b_dummy) {
      com = dummy_atom_pos;
      if (cvm::debug()) {
        cvm::log("Dummy atom center of mass = "+cvm::to_str(com)+"\n");
      }
    } else if (is_enabled(f_ag_scalable)) {
      com = (cvm::proxy)->get_atom_group_com(index);
    } else {
      com.reset();
      for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
     com += ai->mass * ai->pos;     com += ai->mass * ai->pos;
   }   }
   com /= this->total_mass;     com /= total_mass;
   return com;   }
    return COLVARS_OK;
  }
  
  
  int cvm::atom_group::calc_dipole(cvm::atom_pos const &com)
  {
    if (b_dummy) {
      cvm::error("Error: trying to compute the dipole of an empty group.\n", INPUT_ERROR);
      return COLVARS_ERROR;
    }
    dip.reset();
    for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
      dip += ai->charge * (ai->pos - com);
    }
    return COLVARS_OK;
 } }
  
  
Line 646
Line 943
 { {
   if (b_dummy) return;   if (b_dummy) return;
  
   for (cvm::atom_iter ai = this->begin();   if (is_enabled(f_ag_scalable)) {
        ai != this->end(); ai++) {     scalar_com_gradient = grad;
     ai->grad = (ai->mass/this->total_mass) * grad;     return;
    }
  
    for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
      ai->grad = (ai->mass/total_mass) * grad;
   }   }
 } }
  
Line 657
Line 958
 { {
   if (b_dummy) return;   if (b_dummy) return;
  
   if ((!b_center) && (!b_rotate)) return; // no fit 
  
   if (cvm::debug())   if (cvm::debug())
     cvm::log("Calculating fit gradients.\n");     cvm::log("Calculating fit gradients.\n");
  
   atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;   atom_group *group_for_fit = fitting_group ? fitting_group : this;
   group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::rvector(0.0, 0.0, 0.0)); 
  
   if (b_center) {   if (b_center) {
     // add the center of geometry contribution to the gradients     // add the center of geometry contribution to the gradients
      cvm::rvector atom_grad;
  
     for (size_t i = 0; i < this->size(); i++) {     for (size_t i = 0; i < this->size(); i++) {
       // need to bring the gradients in original frame first       atom_grad += atoms[i].grad;
       cvm::rvector const atom_grad = b_rotate ? 
         (rot.inverse()).rotate((*this)[i].grad) : 
         (*this)[i].grad; 
       for (size_t j = 0; j < group_for_fit->size(); j++) { 
         group_for_fit->fit_gradients[j] += 
           (-1.0)/(cvm::real(group_for_fit->size())) * 
           atom_grad; 
       }       }
      if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad);
      atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
  
      for (size_t j = 0; j < group_for_fit->size(); j++) {
        group_for_fit->fit_gradients[j] = atom_grad;
     }     }
   }   }
  
Line 684
Line 982
  
     // add the rotation matrix contribution to the gradients     // add the rotation matrix contribution to the gradients
     cvm::rotation const rot_inv = rot.inverse();     cvm::rotation const rot_inv = rot.inverse();
     cvm::atom_pos const cog = this->center_of_geometry(); 
  
     for (size_t i = 0; i < this->size(); i++) {     for (size_t i = 0; i < this->size(); i++) {
  
       cvm::atom_pos const pos_orig = rot_inv.rotate((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));       // compute centered, unrotated position
        cvm::atom_pos const pos_orig =
          rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
  
       for (size_t j = 0; j < group_for_fit->size(); j++) { 
         // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i         // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
         cvm::quaternion const dxdq =         cvm::quaternion const dxdq =
           rot.q.position_derivative_inner(pos_orig, (*this)[i].grad);         rot.q.position_derivative_inner(pos_orig, atoms[i].grad);
         // multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients 
        for (size_t j = 0; j < group_for_fit->size(); j++) {
          // multiply by {\partial q}/\partial\vec{x}_j and add it to the fit gradients
         for (size_t iq = 0; iq < 4; iq++) {         for (size_t iq = 0; iq < 4; iq++) {
           group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];           group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
         }         }
       }       }
     }     }
   }   }
  
   if (cvm::debug())   if (cvm::debug())
     cvm::log("Done calculating fit gradients.\n");     cvm::log("Done calculating fit gradients.\n");
 } }
Line 710
Line 1011
 { {
   if (b_dummy) {   if (b_dummy) {
     cvm::error("Error: positions are not available "     cvm::error("Error: positions are not available "
                       "from a dummy atom group.\n");                "from a dummy atom group.\n", INPUT_ERROR);
    }
  
    if (is_enabled(f_ag_scalable)) {
      cvm::error("Error: atomic positions are not available "
                 "from a scalable atom group.\n", INPUT_ERROR);
   }   }
  
   std::vector<cvm::atom_pos> x(this->size(), 0.0);   std::vector<cvm::atom_pos> x(this->size(), 0.0);
Line 726
Line 1032
 { {
   if (b_dummy) {   if (b_dummy) {
     cvm::error("Error: positions are not available "     cvm::error("Error: positions are not available "
                  "from a dummy atom group.\n");                "from a dummy atom group.\n", INPUT_ERROR);
    }
  
    if (is_enabled(f_ag_scalable)) {
      cvm::error("Error: atomic positions are not available "
                 "from a scalable atom group.\n", INPUT_ERROR);
   }   }
  
   std::vector<cvm::atom_pos> x(this->size(), 0.0);   std::vector<cvm::atom_pos> x(this->size(), 0.0);
Line 742
Line 1053
 { {
   if (b_dummy) {   if (b_dummy) {
     cvm::error("Error: velocities are not available "     cvm::error("Error: velocities are not available "
                  "from a dummy atom group.\n");                "from a dummy atom group.\n", INPUT_ERROR);
    }
  
    if (is_enabled(f_ag_scalable)) {
      cvm::error("Error: atomic velocities are not available "
                 "from a scalable atom group.\n", INPUT_ERROR);
   }   }
  
   std::vector<cvm::rvector> v(this->size(), 0.0);   std::vector<cvm::rvector> v(this->size(), 0.0);
Line 754
Line 1070
   return v;   return v;
 } }
  
 std::vector<cvm::rvector> cvm::atom_group::system_forces() const std::vector<cvm::rvector> cvm::atom_group::total_forces() const
 { {
   if (b_dummy) {   if (b_dummy) {
     cvm::error("Error: system forces are not available "     cvm::error("Error: total forces are not available "
                  "from a dummy atom group.\n");                "from a dummy atom group.\n", INPUT_ERROR);
    }
  
    if (is_enabled(f_ag_scalable)) {
      cvm::error("Error: atomic total forces are not available "
                 "from a scalable atom group.\n", INPUT_ERROR);
   }   }
  
   std::vector<cvm::rvector> f(this->size(), 0.0);   std::vector<cvm::rvector> f(this->size(), 0.0);
   cvm::atom_const_iter ai = this->begin();   cvm::atom_const_iter ai = this->begin();
   std::vector<cvm::atom_pos>::iterator fi = f.begin();   std::vector<cvm::atom_pos>::iterator fi = f.begin();
   for ( ; ai != this->end(); ++fi, ++ai) {   for ( ; ai != this->end(); ++fi, ++ai) {
     *fi = ai->system_force;     *fi = ai->total_force;
   }   }
   return f;   return f;
 } }
  
 cvm::rvector cvm::atom_group::system_force() const 
  // TODO make this an accessor
  cvm::rvector cvm::atom_group::total_force() const
 { {
   if (b_dummy) {   if (b_dummy) {
     cvm::error("Error: system forces are not available "     cvm::error("Error: total total forces are not available "
                 "from a dummy atom group.\n");                "from a dummy atom group.\n", INPUT_ERROR);
    }
  
    if (is_enabled(f_ag_scalable)) {
      return (cvm::proxy)->get_atom_group_total_force(index);
   }   }
  
   cvm::rvector f(0.0);   cvm::rvector f(0.0);
   for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {   for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
     f += ai->system_force;     f += ai->total_force;
   }   }
   return f;   return f;
 } }
Line 787
Line 1114
  
 void cvm::atom_group::apply_colvar_force(cvm::real const &force) void cvm::atom_group::apply_colvar_force(cvm::real const &force)
 { {
    if (cvm::debug()) {
      log("Communicating a colvar force from atom group to the MD engine.\n");
    }
  
   if (b_dummy)   if (b_dummy)
     return;     return;
  
Line 796
Line 1127
     return;     return;
   }   }
  
    if (is_enabled(f_ag_scalable)) {
      (cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient);
      return;
    }
  
   if (b_rotate) {   if (b_rotate) {
  
     // rotate forces back to the original frame     // rotate forces back to the original frame
     cvm::rotation const rot_inv = rot.inverse();     cvm::rotation const rot_inv = rot.inverse();
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) { 
       ai->apply_force(rot_inv.rotate(force * ai->grad));       ai->apply_force(rot_inv.rotate(force * ai->grad));
     }     }
  
   } else {   } else {
  
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) { 
       ai->apply_force(force * ai->grad);       ai->apply_force(force * ai->grad);
     }     }
   }   }
  
   if ((b_center || b_rotate) && b_fit_gradients) {   if ((b_center || b_rotate) && b_fit_gradients) {
  
     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;     atom_group *group_for_fit = fitting_group ? fitting_group : this;
  
     // add the contribution from the roto-translational fit to the gradients     // Fit gradients are already calculated in "laboratory" frame
     if (b_rotate) { 
       // rotate forces back to the original frame 
       cvm::rotation const rot_inv = rot.inverse(); 
       for (size_t j = 0; j < group_for_fit->size(); j++) { 
         (*group_for_fit)[j].apply_force(rot_inv.rotate(force * group_for_fit->fit_gradients[j])); 
       } 
     } else { 
       for (size_t j = 0; j < group_for_fit->size(); j++) {       for (size_t j = 0; j < group_for_fit->size(); j++) {
         (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);         (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
       }       }
     }     }
   }   }
  
 } 
  
  
 void cvm::atom_group::apply_force(cvm::rvector const &force) void cvm::atom_group::apply_force(cvm::rvector const &force)
 { {
Line 845
Line 1170
     return;     return;
   }   }
  
    if (is_enabled(f_ag_scalable)) {
      (cvm::proxy)->apply_atom_group_force(index, force);
    }
  
   if (b_rotate) {   if (b_rotate) {
  
     cvm::rotation const rot_inv = rot.inverse();     cvm::rotation const rot_inv = rot.inverse();
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) {       ai->apply_force(rot_inv.rotate((ai->mass/total_mass) * force));
       ai->apply_force(rot_inv.rotate((ai->mass/this->total_mass) * force)); 
     }     }
  
   } else {   } else {
  
     for (cvm::atom_iter ai = this->begin();     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
          ai != this->end(); ai++) {       ai->apply_force((ai->mass/total_mass) * force);
       ai->apply_force((ai->mass/this->total_mass) * force); 
     }     }
   }   }
 } }
  
  
 void cvm::atom_group::apply_forces(std::vector<cvm::rvector> const &forces) // Static members
 { 
   if (b_dummy) 
     return; 
  
   if (noforce) std::vector<colvardeps::feature *> cvm::atom_group::ag_features;
     cvm::error("Error: sending a force to a group that has " 
                 "\"disableForces\" defined.\n"); 
  
   if (forces.size() != this->size()) { 
     cvm::error("Error: trying to apply an array of forces to an atom " 
                 "group which does not have the same length.\n"); 
   } 
  
   if (b_rotate) { 
  
     cvm::rotation const rot_inv = rot.inverse(); 
     cvm::atom_iter ai = this->begin(); 
     std::vector<cvm::rvector>::const_iterator fi = forces.begin(); 
     for ( ; ai != this->end(); ++fi, ++ai) { 
       ai->apply_force(rot_inv.rotate(*fi)); 
     } 
  
   } else { 
  
     cvm::atom_iter ai = this->begin(); 
     std::vector<cvm::rvector>::const_iterator fi = forces.begin(); 
     for ( ; ai != this->end(); ++fi, ++ai) { 
       ai->apply_force(*fi); 
     } 
   } 
 } 
  


Legend:
Removed in v.1.19 
changed lines
 Added in v.1.20



Made by using version 1.53 of cvs2html