Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision d103d335f1f411ee26c2568846556af546c71b43)
+++ src/molecule.cpp	(revision 30c7536948400224f2fc3e945a95099d12c23b3a)
@@ -150,22 +150,22 @@
 
 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());
 }
 
@@ -182,8 +182,7 @@
   OBSERVE;
   molecule::const_iterator iter = loc;
-  iter++;
+  ++iter;
   atom* atom = *loc;
   atomIds.erase( atom->getId() );
-  atoms.remove( atom );
   formula-=atom->getType();
   atom->removeFromMolecule();
@@ -196,7 +195,6 @@
   molecule::const_iterator iter = find(key);
   if (iter != end()){
-    iter++;
+    ++iter;
     atomIds.erase( key->getId() );
-    atoms.remove( key );
     formula-=key->getType();
     key->removeFromMolecule();
@@ -207,10 +205,5 @@
 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,9 +213,8 @@
   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,6 +227,12 @@
 {
   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::const_iterator iter = begin(); iter != end(); ++iter)
+    vector_of_atoms.push_back(*iter);
   return vector_of_atoms;
 }
@@ -635,5 +633,67 @@
 
   // copy all atoms
-  for_each(atoms.begin(),atoms.end(),bind1st(mem_fun(&molecule::AddCopyAtom),copy));
+  std::map< const atom *, atom *> FatherFinder;
+  for (const_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 for bond copy found");
+        ASSERT(FatherFinder.count(Binder->rightatom),
+            "molecule::CopyMolecule() - No copy of original right atom 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(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();
+
+  // copy all atoms
+  std::map< const atom *, atom *> FatherFinder;
+  for (const_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,50 +705,20 @@
       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);
@@ -878,5 +908,5 @@
     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,5 +930,5 @@
       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,5 +971,5 @@
     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,5 +988,5 @@
     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,5 +1002,5 @@
   int i = 0;
   NoNonHydrogen = 0;
-  for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+  for (molecule::const_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
