Changes in src/molecule.cpp [d103d3:59fff1]

File:

• src/molecule.cpp (modified) (16 diffs)

src/molecule.cpp

@@ -150 +150 @@
 
 molecule::iterator molecule::begin(){
-  return molecule::iterator(atoms.begin(),this);
+  return iterator(atomIds.begin(), FromIdToAtom());
 }
 
 molecule::const_iterator molecule::begin() const{
-  return atoms.begin();
+  return const_iterator(atomIds.begin(), FromIdToAtom());
 }
 
 molecule::iterator molecule::end(){
-  return molecule::iterator(atoms.end(),this);
+  return iterator(atomIds.end(), FromIdToAtom());
 }
 
 molecule::const_iterator molecule::end() const{
-  return atoms.end();
+  return const_iterator(atomIds.end(), FromIdToAtom());
 }
 
 bool molecule::empty() const
 {
-  return (begin() == end());
+  return (atomIds.empty());
 }
 
@@ -173 +173 @@
 {
   size_t counter = 0;
-  for (molecule::const_iterator iter = begin(); iter != end (); ++iter)
+  for (const_iterator iter = begin(); iter != end (); ++iter)
     counter++;
   return counter;
@@ -181 +181 @@
 {
   OBSERVE;
-  molecule::const_iterator iter = loc;
-  iter++;
-  atom* atom = *loc;
-  atomIds.erase( atom->getId() );
-  atoms.remove( atom );
-  formula-=atom->getType();
-  atom->removeFromMolecule();
+  const_iterator iter = loc;
+  ++iter;
+  atom * const _atom = const_cast<atom *>(*loc);
+  atomIds.erase( _atom->getId() );
+  formula-=_atom->getType();
+  _atom->removeFromMolecule();
   return iter;
 }
@@ -194 +193 @@
 {
   OBSERVE;
-  molecule::const_iterator iter = find(key);
+  const_iterator iter = find(key);
   if (iter != end()){
-    iter++;
+    ++iter;
     atomIds.erase( key->getId() );
-    atoms.remove( key );
     formula-=key->getType();
     key->removeFromMolecule();
@@ -207 +205 @@
 molecule::const_iterator molecule::find ( atom * key ) const
 {
-  molecule::const_iterator iter;
-  for (molecule::const_iterator Runner = begin(); Runner != end(); ++Runner) {
-    if (*Runner == key)
-      return molecule::const_iterator(Runner);
-  }
-  return molecule::const_iterator(atoms.end());
+  return const_iterator(atomIds.find(key->getId()), FromIdToAtom());
 }
 
@@ -220 +213 @@
   pair<atomIdSet::iterator,bool> res = atomIds.insert(key->getId());
   if (res.second) { // push atom if went well
-    atoms.push_back(key);
     formula+=key->getType();
-    return pair<iterator,bool>(molecule::iterator(--end()),res.second);
+    return pair<iterator,bool>(iterator(res.first, FromIdToAtom()),res.second);
   } else {
-    return pair<iterator,bool>(molecule::iterator(end()),res.second);
+    return pair<iterator,bool>(end(),res.second);
   }
 }
@@ -235 +227 @@
 {
   World::AtomComposite vector_of_atoms;
-  BOOST_FOREACH(atom *_atom, atoms)
-    vector_of_atoms.push_back(_atom);
+//  std::copy(MyIter(atomIds.begin(), FromIdToAtom()),
+//      MyIter(atomIds.end(), FromIdToAtom()),
+//      vector_of_atoms.begin());
+//  for (MyIter iter = MyIter(atomIds.begin(), FromIdToAtom());
+//      iter != MyIter(atomIds.end(), FromIdToAtom());
+//      ++iter)
+  for (molecule::iterator iter = begin(); iter != end(); ++iter)
+    vector_of_atoms.push_back(*iter);
   return vector_of_atoms;
 }
@@ -635 +633 @@
 
   // copy all atoms
-  for_each(atoms.begin(),atoms.end(),bind1st(mem_fun(&molecule::AddCopyAtom),copy));
+  std::map< const atom *, atom *> FatherFinder;
+  for (iterator iter = begin(); iter != end(); ++iter) {
+    atom * const copy_atom = copy->AddCopyAtom(*iter);
+    FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
+  }
+
+  // copy all bonds
+  for(const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
+      if ((*BondRunner)->leftatom == *AtomRunner) {
+        bond *Binder = (*BondRunner);
+        // get the pendant atoms of current bond in the copy molecule
+        ASSERT(FatherFinder.count(Binder->leftatom),
+            "molecule::CopyMolecule() - No copy of original left atom "
+            +toString(Binder->leftatom)+" for bond copy found");
+        ASSERT(FatherFinder.count(Binder->rightatom),
+            "molecule::CopyMolecule() - No copy of original right atom "
+            +toString(Binder->rightatom)+"  for bond copy found");
+        atom * const LeftAtom = FatherFinder[Binder->leftatom];
+        atom * const RightAtom = FatherFinder[Binder->rightatom];
+
+        bond * const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
+        NewBond->Cyclic = Binder->Cyclic;
+        if (Binder->Cyclic)
+          copy->NoCyclicBonds++;
+        NewBond->Type = Binder->Type;
+      }
+  }
+  // correct fathers
+  //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
+
+  return copy;
+};
+
+
+/** Destroys all atoms inside this molecule.
+ */
+void molecule::removeAtomsinMolecule()
+{
+  // remove each atom from world
+  for(iterator AtomRunner = begin(); !empty(); AtomRunner = begin())
+    World::getInstance().destroyAtom(*AtomRunner);
+};
+
+
+/**
+ * Copies all atoms of a molecule which are within the defined parallelepiped.
+ *
+ * @param offest for the origin of the parallelepiped
+ * @param three vectors forming the matrix that defines the shape of the parallelpiped
+ */
+molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
+  molecule *copy = World::getInstance().createMolecule();
+
+  // copy all atoms
+  std::map< const atom *, atom *> FatherFinder;
+  for (iterator iter = begin(); iter != end(); ++iter) {
+    if((*iter)->IsInShape(region)){
+      atom * const copy_atom = copy->AddCopyAtom(*iter);
+      FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
+    }
+  }
 
   // copy all bonds
@@ -645 +707 @@
       if ((*BondRunner)->leftatom == *AtomRunner) {
         bond *Binder = (*BondRunner);
-        // get the pendant atoms of current bond in the copy molecule
-        atomSet::iterator leftiter=find_if(copy->atoms.begin(),copy->atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->leftatom));
-        atomSet::iterator rightiter=find_if(copy->atoms.begin(),copy->atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->rightatom));
-        ASSERT(leftiter!=copy->atoms.end(),"No copy of original left atom for bond copy found");
-        ASSERT(leftiter!=copy->atoms.end(),"No copy of original right atom for bond copy found");
-        atom *LeftAtom = *leftiter;
-        atom *RightAtom = *rightiter;
-
-        bond *NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
-        NewBond->Cyclic = Binder->Cyclic;
-        if (Binder->Cyclic)
-          copy->NoCyclicBonds++;
-        NewBond->Type = Binder->Type;
+        if ((FatherFinder.count(Binder->leftatom))
+            && (FatherFinder.count(Binder->rightatom))) {
+          // if copy present, then it must be from subregion
+          atom * const LeftAtom = FatherFinder[Binder->leftatom];
+          atom * const RightAtom = FatherFinder[Binder->rightatom];
+
+          bond * const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
+          NewBond->Cyclic = Binder->Cyclic;
+          if (Binder->Cyclic)
+            copy->NoCyclicBonds++;
+          NewBond->Type = Binder->Type;
+        }
       }
   }
   // correct fathers
-  //for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
-
-  return copy;
-};
-
-
-/** Destroys all atoms inside this molecule.
- */
-void molecule::removeAtomsinMolecule()
-{
-  // remove each atom from world
-  for(molecule::const_iterator AtomRunner = begin(); !empty(); AtomRunner = begin())
-    World::getInstance().destroyAtom(*AtomRunner);
-};
-
-
-/**
- * Copies all atoms of a molecule which are within the defined parallelepiped.
- *
- * @param offest for the origin of the parallelepiped
- * @param three vectors forming the matrix that defines the shape of the parallelpiped
- */
-molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
-  molecule *copy = World::getInstance().createMolecule();
-
-  BOOST_FOREACH(atom *iter,atoms){
-    if(iter->IsInShape(region)){
-      copy->AddCopyAtom(iter);
-    }
-  }
+  //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
 
   //TODO: copy->BuildInducedSubgraph(this);
@@ -709 +741 @@
 
   // some checks to make sure we are able to create the bond
-  ASSERT(atom1, "First atom in bond-creation was an invalid pointer");
-  ASSERT(atom2, "Second atom in bond-creation was an invalid pointer");
-  ASSERT(FindAtom(atom1->getNr()),"First atom in bond-creation was not part of molecule");
-  ASSERT(FindAtom(atom2->getNr()),"Second atom in bond-creation was not part of molecule");
+  ASSERT(atom1,
+      "molecule::AddBond() - First atom "+toString(atom1)
+      +" is not a invalid pointer");
+  ASSERT(atom2,
+      "molecule::AddBond() - Second atom "+toString(atom2)
+      +" is not a invalid pointer");
+  ASSERT(isInMolecule(atom1),
+      "molecule::AddBond() - First atom "+toString(atom1)
+      +" is not part of molecule");
+  ASSERT(isInMolecule(atom2),
+      "molecule::AddBond() - Second atom "+toString(atom2)
+      +" is not part of molecule");
 
   Binder = new bond(atom1, atom2, degree);
   atom1->RegisterBond(WorldTime::getTime(), Binder);
   atom2->RegisterBond(WorldTime::getTime(), Binder);
-  if ((atom1->getType() != NULL) && (atom1->getType()->getAtomicNumber() != 1) && (atom2->getType() != NULL) && (atom2->getType()->getAtomicNumber() != 1))
+  if ((atom1->getType() != NULL)
+      && (atom1->getType()->getAtomicNumber() != 1)
+      && (atom2->getType() != NULL)
+      && (atom2->getType()->getAtomicNumber() != 1))
     NoNonBonds++;
 
@@ -820 +863 @@
 atom * molecule::FindAtom(int Nr)  const
 {
-  molecule::const_iterator iter = begin();
+  molecule::iterator iter = begin();
   for (; iter != end(); ++iter)
-    if ((*iter)->getNr() == Nr)
-      break;
+  if ((*iter)->getNr() == Nr)
+    break;
   if (iter != end()) {
     //LOG(0, "Found Atom Nr. " << walker->getNr());
     return (*iter);
   } else {
-    LOG(0, "Atom not found in list.");
+    ELOG(1, "Atom not found in molecule " << getName() << "'s list.");
     return NULL;
   }
-};
+}
+
+/** Checks whether the given atom is a member of this molecule.
+ *
+ *  We make use here of molecule::atomIds to get a result on
+ *
+ * @param _atom atom to check
+ * @return true - is member, false - is not
+ */
+bool molecule::isInMolecule(const atom * const _atom)
+{
+  ASSERT(_atom->getMolecule() == this,
+      "molecule::isInMolecule() - atom is not designated to be in molecule '"
+      +toString(this->getName())+"'.");
+  molecule::atomIdSet::const_iterator iter = atomIds.find(_atom->getId());
+  return (iter != atomIds.end());
+}
 
 /** Asks for atom number, and checks whether in list.
@@ -878 +937 @@
     enumeration<const element*> elementLookup = formula.enumerateElements();
     *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
-    for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
+    for_each(begin(),end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
     return true;
   }
@@ -900 +959 @@
       memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
       enumeration<const element*> elementLookup = formula.enumerateElements();
-      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectory,_1,output,elementLookup, AtomNo, (const int)step));
+      for_each(begin(),end(),boost::bind(&atom::OutputTrajectory,_1,output,elementLookup, AtomNo, (const int)step));
     }
     return true;
@@ -941 +1000 @@
     for (int step=0;step<MDSteps;step++) {
       *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
-      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
+      for_each(begin(),end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
     }
     return true;
@@ -958 +1017 @@
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
-    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
+    for_each(begin(),end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
     return true;
   } else
@@ -972 +1031 @@
   int i = 0;
   NoNonHydrogen = 0;
-  for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+  for (molecule::iterator iter = begin(); iter != end(); ++iter) {
     (*iter)->setNr(i);   // update number in molecule (for easier referencing in FragmentMolecule lateron)
     if ((*iter)->getType()->getAtomicNumber() != 1) // count non-hydrogen atoms whilst at it