Index: src/Legacy/oldmenu.cpp
===================================================================
--- src/Legacy/oldmenu.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/Legacy/oldmenu.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -447,14 +447,14 @@
       break;
     case 'b':
-      second = mol->AskAtom("Enter number of atom as reference point: ");
-      Log() << Verbose(0) << "Enter radius: ";
-      cin >> tmp1;
-      first = mol->start;
-      second = first->next;
-      while(second != mol->end) {
-        first = second;
-        second = first->next;
-        if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
-          mol->RemoveAtom(first);
+      {
+        second = mol->AskAtom("Enter number of atom as reference point: ");
+        Log() << Verbose(0) << "Enter radius: ";
+        cin >> tmp1;
+        molecule::iterator runner;
+        for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
+          runner = iter++;
+          if ((*runner)->x.DistanceSquared((const Vector *)&(*runner)->x) > tmp1*tmp1) // distance to first above radius ...
+            mol->RemoveAtom((*runner));
+        }
       }
       break;
@@ -466,12 +466,10 @@
       Log() << Verbose(0) << "Upper boundary: ";
       cin >> tmp2;
-      first = mol->start;
-      second = first->next;
-      while(second != mol->end) {
-        first = second;
-        second = first->next;
-        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
-          //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
-          mol->RemoveAtom(first);
+      molecule::iterator runner;
+      for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
+        runner = iter++;
+        if (((*runner)->x.x[axis] < tmp1) || ((*runner)->x.x[axis] > tmp2)) {// out of boundary ...
+          //Log() << Verbose(0) << "Atom " << *(*runner) << " with " << (*runner)->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
+          mol->RemoveAtom((*runner));
         }
       }
@@ -516,16 +514,14 @@
         min[i] = 0.;
 
-      second = mol->start;
-      while ((second->next != mol->end)) {
-        second = second->next; // advance
-        Z = second->type->Z;
+      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+        Z = (*iter)->type->Z;
         tmp1 = 0.;
-        if (first != second) {
+        if (first != (*iter)) {
           x.CopyVector((const Vector *)&first->x);
-          x.SubtractVector((const Vector *)&second->x);
+          x.SubtractVector((const Vector *)&(*iter)->x);
           tmp1 = x.Norm();
         }
         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
-        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
+        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << ((*iter)->nr << ": " << tmp1 << " a.u." << endl;
       }
       for (int i=MAX_ELEMENTS;i--;)
@@ -767,9 +763,7 @@
       vectors = new Vector *[count];
       j = 0;
-      first = mol->start;
-      while (first->next != mol->end) { // make a list of all atoms with coordinates and element
-        first = first->next;
-        Elements[j] = first->type;
-        vectors[j] = &first->x;
+      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+        Elements[j] = (*iter)->type;
+        vectors[j] = &(*iter)->x;
         j++;
       }
@@ -1029,5 +1023,4 @@
     return;
   }
-  atom *Walker = mol->start;
 
   // generate some KeySets
@@ -1035,8 +1028,7 @@
   KeySet TestSets[mol->AtomCount+1];
   i=1;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     for (int j=0;j<i;j++) {
-      TestSets[j].insert(Walker->nr);
+      TestSets[j].insert((*iter)->nr);
     }
     i++;
@@ -1044,12 +1036,15 @@
   Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
   KeySetTestPair test;
-  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
-  if (test.second) {
-    Log() << Verbose(1) << "Insertion worked?!" << endl;
+  molecule::const_iterator iter = mol->begin();
+  if (iter != mol->end()) {
+    test = TestSets[mol->AtomCount-1].insert((*iter)->nr);
+    if (test.second) {
+      Log() << Verbose(1) << "Insertion worked?!" << endl;
+    } else {
+      Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
+    }
   } else {
-    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
-  }
-  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
-  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
+    eLog() << Verbose(1) << "No atoms to test double insertion." << endl;
+  }
 
   // constructing Graph structure
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/Makefile.am	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -41,5 +41,5 @@
 				   Descriptors/MoleculeIdDescriptor.hpp
 
-SOURCE = ${ANALYSISSOURCE} ${ATOMSOURCE} ${PATTERNSOURCE} ${UISOURCE} ${DESCRIPTORSOURCE} ${LEGACYSOURCE} bond.cpp bondgraph.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp errorlogger.cpp graph.cpp helpers.cpp info.cpp leastsquaremin.cpp linkedcell.cpp lists.cpp log.cpp logger.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp World.cpp WorldIterators.cpp
+SOURCE = ${ANALYSISSOURCE} ${ATOMSOURCE} ${PATTERNSOURCE} ${UISOURCE} ${DESCRIPTORSOURCE} ${LEGACYSOURCE} bond.cpp bondgraph.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp errorlogger.cpp graph.cpp helpers.cpp info.cpp leastsquaremin.cpp linkedcell.cpp log.cpp logger.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp World.cpp WorldIterators.cpp
 HEADER = ${ANALYSISHEADER} ${ATOMHEADER} ${PATTERNHEADER} ${UIHEADER} ${DESCRIPTORHEADER} ${LEGACYHEADER} bond.hpp bondgraph.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp errorlogger.hpp graph.hpp helpers.hpp info.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp log.hpp logger.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp World.hpp 
 
Index: src/analysis_bonds.cpp
===================================================================
--- src/analysis_bonds.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/analysis_bonds.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -24,9 +24,7 @@
   Mean = 0.;
 
-  atom *Walker = mol->start;
   int AtomCount = 0;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    const int count = Walker->ListOfBonds.size();
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    const int count = (*iter)->ListOfBonds.size();
     if (Max < count)
       Max = count;
@@ -56,10 +54,8 @@
 
   int AtomNo = 0;
-  atom *Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    if (Walker->type == type1)
-      for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++)
-        if ((*BondRunner)->GetOtherAtom(Walker)->type == type2) {
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    if ((*iter)->type == type1)
+      for (BondList::const_iterator BondRunner = (*iter)->ListOfBonds.begin(); BondRunner != (*iter)->ListOfBonds.end(); BondRunner++)
+        if ((*BondRunner)->GetOtherAtom((*iter))->type == type2) {
           const double distance = (*BondRunner)->GetDistanceSquared();
           if (Min > distance)
Index: src/analysis_correlation.cpp
===================================================================
--- src/analysis_correlation.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/analysis_correlation.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -41,21 +41,17 @@
     if ((*MolWalker)->ActiveFlag) {
       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
-        if ((type1 == NULL) || (Walker->type == type1)) {
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
+        if ((type1 == NULL) || ((*iter)->type == type1)) {
           for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
             if ((*MolOtherWalker)->ActiveFlag) {
               Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
-              atom *OtherWalker = (*MolOtherWalker)->start;
-              while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
-                OtherWalker = OtherWalker->next;
-                Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
-                if (Walker->nr < OtherWalker->nr)
-                  if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
-                    //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
-                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
+              for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
+                Log() << Verbose(3) << "Current otheratom is " << *(*runner) << "." << endl;
+                if ((*iter)->nr < (*runner)->nr)
+                  if ((type2 == NULL) || ((*runner)->type == type2)) {
+                    distance = (*iter)->node->PeriodicDistance((*runner)->node, (*MolWalker)->cell_size);
+                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
+                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
                   }
               }
@@ -99,10 +95,8 @@
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
-        if ((type1 == NULL) || (Walker->type == type1)) {
-          periodicX.CopyVector(Walker->node);
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
+        if ((type1 == NULL) || ((*iter)->type == type1)) {
+          periodicX.CopyVector((*iter)->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
@@ -116,11 +110,9 @@
                   if ((*MolOtherWalker)->ActiveFlag) {
                     Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
-                    atom *OtherWalker = (*MolOtherWalker)->start;
-                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
-                      OtherWalker = OtherWalker->next;
-                      Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
-                      if (Walker->nr < OtherWalker->nr)
-                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                          periodicOtherX.CopyVector(OtherWalker->node);
+                    for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
+                      Log() << Verbose(3) << "Current otheratom is " << *(*runner) << "." << endl;
+                      if ((*iter)->nr < (*runner)->nr)
+                        if ((type2 == NULL) || ((*runner)->type == type2)) {
+                          periodicOtherX.CopyVector((*runner)->node);
                           periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
                           // go through every range in xyz and get distance
@@ -132,6 +124,6 @@
                                 checkOtherX.MatrixMultiplication(FullMatrix);
                                 distance = checkX.Distance(&checkOtherX);
-                                //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
-                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
+                                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
+                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
                               }
                         }
@@ -169,12 +161,10 @@
     if ((*MolWalker)->ActiveFlag) {
       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
-        if ((type == NULL) || (Walker->type == type)) {
-          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
+        if ((type == NULL) || ((*iter)->type == type)) {
+          distance = (*iter)->node->PeriodicDistance(point, (*MolWalker)->cell_size);
           Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
-          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
+          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
         }
       }
@@ -211,10 +201,8 @@
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
-        if ((type == NULL) || (Walker->type == type)) {
-          periodicX.CopyVector(Walker->node);
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
+        if ((type == NULL) || ((*iter)->type == type)) {
+          periodicX.CopyVector((*iter)->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
@@ -227,5 +215,5 @@
                 distance = checkX.Distance(point);
                 Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
-                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
+                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
               }
         }
@@ -262,13 +250,11 @@
     if ((*MolWalker)->ActiveFlag) {
       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
-        if ((type == NULL) || (Walker->type == type)) {
-          triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
+        if ((type == NULL) || ((*iter)->type == type)) {
+          triangle = Surface->FindClosestTriangleToVector((*iter)->node, LC );
           if (triangle != NULL) {
-            distance = DistanceToTrianglePlane(Walker->node, triangle);
-            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
+            distance = DistanceToTrianglePlane((*iter)->node, triangle);
+            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
           }
         }
@@ -315,10 +301,8 @@
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
-        if ((type == NULL) || (Walker->type == type)) {
-          periodicX.CopyVector(Walker->node);
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
+        if ((type == NULL) || ((*iter)->type == type)) {
+          periodicX.CopyVector((*iter)->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
@@ -339,6 +323,6 @@
           // insert
           ShortestDistance = sqrt(ShortestDistance);
-          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
-          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
+          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> ((*iter), ShortestTriangle) ) );
+          //Log() << Verbose(1) << "INFO: Inserting " << (*iter) << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
         }
       }
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/atom.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -59,5 +59,4 @@
 atom::~atom()
 {
-  unlink(this);
 };
 
Index: src/bondgraph.cpp
===================================================================
--- src/bondgraph.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/bondgraph.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -88,5 +88,5 @@
   bool status = true;
 
-  if (mol->start->next == mol->end) // only construct if molecule is not empty
+  if (mol->empty()) // only construct if molecule is not empty
     return false;
 
@@ -121,9 +121,7 @@
   max_distance = 0.;
 
-  atom *Runner = mol->start;
-  while (Runner->next != mol->end) {
-    Runner = Runner->next;
-    if (Runner->type->CovalentRadius > max_distance)
-      max_distance = Runner->type->CovalentRadius;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    if ((*iter)->type->CovalentRadius > max_distance)
+      max_distance = (*iter)->type->CovalentRadius;
   }
   max_distance *= 2.;
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/boundary.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -140,5 +140,4 @@
 {
 	Info FunctionInfo(__func__);
-  atom *Walker = NULL;
   PointMap PointsOnBoundary;
   LineMap LinesOnBoundary;
@@ -166,8 +165,6 @@
 
     // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
-    Walker = mol->start;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      ProjectedVector.CopyVector(&Walker->x);
+    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      ProjectedVector.CopyVector(&(*iter)->x);
       ProjectedVector.SubtractVector(MolCenter);
       ProjectedVector.ProjectOntoPlane(&AxisVector);
@@ -184,17 +181,17 @@
         angle = 2. * M_PI - angle;
       }
-      Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
-      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
+      Log() << Verbose(1) << "Inserting " << *(*iter) << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
+      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, (*iter))));
       if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
         Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
         Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
-        Log() << Verbose(2) << "New vector: " << *Walker << endl;
+        Log() << Verbose(2) << "New vector: " << *(*iter) << endl;
         const double ProjectedVectorNorm = ProjectedVector.NormSquared();
         if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
           BoundaryTestPair.first->second.first = ProjectedVectorNorm;
-          BoundaryTestPair.first->second.second = Walker;
+          BoundaryTestPair.first->second.second = (*iter);
           Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
         } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper.CopyVector(&Walker->x);
+          helper.CopyVector(&(*iter)->x);
           helper.SubtractVector(MolCenter);
           const double oldhelperNorm = helper.NormSquared();
@@ -202,5 +199,5 @@
           helper.SubtractVector(MolCenter);
           if (helper.NormSquared() < oldhelperNorm) {
-            BoundaryTestPair.first->second.second = Walker;
+            BoundaryTestPair.first->second.second = (*iter);
             Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
           } else {
@@ -705,5 +702,4 @@
   int repetition[NDIM] = { 1, 1, 1 };
   int TotalNoClusters = 1;
-  atom *Walker = NULL;
   double totalmass = 0.;
   double clustervolume = 0.;
@@ -729,8 +725,6 @@
 
   // sum up the atomic masses
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      totalmass += Walker->type->mass;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      totalmass += (*iter)->type->mass;
   }
   Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
@@ -811,5 +805,4 @@
   Vector FillerDistance;
   double FillIt = false;
-  atom *Walker = NULL;
   bond *Binder = NULL;
   int i = 0;
@@ -884,9 +877,7 @@
 
           // go through all atoms
-          Walker = filler->start;
-          while (Walker->next != filler->end) {
-            Walker = Walker->next;
+          for (molecule::const_iterator iter = filler->begin(); iter != filler->end(); ++iter) {
             // copy atom ...
-            CopyAtoms[Walker->nr] = Walker->clone();
+            CopyAtoms[(*iter)->nr] = (*iter)->clone();
 
             // create atomic random translation vector ...
@@ -913,14 +904,14 @@
             // ... and put at new position
             if (DoRandomRotation)
-              CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
-            CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
-            CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
-            CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
+              CopyAtoms[(*iter)->nr]->x.MatrixMultiplication(Rotations);
+            CopyAtoms[(*iter)->nr]->x.AddVector(&AtomTranslations);
+            CopyAtoms[(*iter)->nr]->x.AddVector(&FillerTranslations);
+            CopyAtoms[(*iter)->nr]->x.AddVector(&CurrentPosition);
 
             // insert into Filling
 
             // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
-            Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
-            Filling->AddAtom(CopyAtoms[Walker->nr]);
+            Log() << Verbose(4) << "Filling atom " << *(*iter) << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->nr]->x) << "." << endl;
+            Filling->AddAtom(CopyAtoms[(*iter)->nr]);
           }
 
@@ -1036,8 +1027,7 @@
 //  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
 //  //->InsertStraddlingPoints(mol, LCList);
-//  mol->GoToFirst();
+//  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
 //  class TesselPoint *Runner = NULL;
-//  while (!mol->IsEnd()) {
-//    Runner = mol->GetPoint();
+//    Runner = *iter;
 //    Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
 //    if (!->IsInnerPoint(Runner, LCList)) {
@@ -1047,5 +1037,4 @@
 //      Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
 //    }
-//    mol->GoToNext();
 //  }
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/builder.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -483,11 +483,9 @@
       Log() << Verbose(0) << "Enter radius: ";
       cin >> tmp1;
-      first = mol->start;
-      second = first->next;
-      while(second != mol->end) {
-        first = second;
-        second = first->next;
-        if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
-          mol->RemoveAtom(first);
+      molecule::iterator runner;
+      for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+        runner = iter++;
+        if ((*iter)->x.DistanceSquared((const Vector *)&(*runner)->x) > tmp1*tmp1) // distance to first above radius ...
+          mol->RemoveAtom((*runner));
       }
       break;
@@ -499,12 +497,10 @@
       Log() << Verbose(0) << "Upper boundary: ";
       cin >> tmp2;
-      first = mol->start;
-      second = first->next;
-      while(second != mol->end) {
-        first = second;
-        second = first->next;
-        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
-          //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
-          mol->RemoveAtom(first);
+      molecule::iterator runner;
+      for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+        runner = iter++;
+        if (((*runner)->x.x[axis] < tmp1) || ((*runner)->x.x[axis] > tmp2)) {// out of boundary ...
+          //Log() << Verbose(0) << "Atom " << *(*runner) << " with " << (*runner)->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
+          mol->RemoveAtom((*runner));
         }
       }
@@ -549,16 +545,14 @@
         min[i] = 0.;
 
-      second = mol->start;
-      while ((second->next != mol->end)) {
-        second = second->next; // advance
-        Z = second->type->Z;
+      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+        Z = (*iter)->type->Z;
         tmp1 = 0.;
-        if (first != second) {
+        if (first != (*iter)) {
           x.CopyVector((const Vector *)&first->x);
-          x.SubtractVector((const Vector *)&second->x);
+          x.SubtractVector((const Vector *)&(*iter)->x);
           tmp1 = x.Norm();
         }
         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
-        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
+        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << (*iter)->nr << ": " << tmp1 << " a.u." << endl;
       }
       for (int i=MAX_ELEMENTS;i--;)
@@ -829,9 +823,7 @@
           vectors = new Vector *[count];
           j = 0;
-          first = mol->start;
-          while (first->next != mol->end) { // make a list of all atoms with coordinates and element
-            first = first->next;
-            Elements[j] = first->type;
-            vectors[j] = &first->x;
+          for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+            Elements[j] = (*iter)->type;
+            vectors[j] = &(*iter)->x;
             j++;
           }
@@ -1184,5 +1176,4 @@
     return;
   }
-  atom *Walker = mol->start;
 
   // generate some KeySets
@@ -1190,8 +1181,7 @@
   KeySet TestSets[mol->AtomCount+1];
   i=1;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     for (int j=0;j<i;j++) {
-      TestSets[j].insert(Walker->nr);
+      TestSets[j].insert((*iter)->nr);
     }
     i++;
@@ -1199,12 +1189,15 @@
   Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
   KeySetTestPair test;
-  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
-  if (test.second) {
-    Log() << Verbose(1) << "Insertion worked?!" << endl;
+  molecule::const_iterator runner = mol->begin();
+  if (runner != mol->end()) {
+    test = TestSets[mol->AtomCount-1].insert((*runner)->nr);
+    if (test.second) {
+      Log() << Verbose(1) << "Insertion worked?!" << endl;
+    } else {
+      Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
+    }
   } else {
-    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
+    eLog() << Verbose(1) << "No atoms to test on." << endl;
   }
-  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
-  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
 
   // constructing Graph structure
@@ -1810,12 +1803,10 @@
                 double tmp1 = atof(argv[argptr+1]);
                 atom *third = mol->FindAtom(atoi(argv[argptr]));
-                atom *first = mol->start;
-                if ((third != NULL) && (first != mol->end)) {
-                  atom *second = first->next;
-                  while(second != mol->end) {
-                    first = second;
-                    second = first->next;
-                    if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
-                      mol->RemoveAtom(first);
+                molecule::iterator runner;
+                if (third != NULL) {
+                  for (molecule::iterator iter = mol->begin(); iter != mol->end();) {
+                    runner = iter++;
+                    if ((*runner)->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
+                      mol->RemoveAtom((*runner));
                   }
                 } else {
@@ -2082,9 +2073,7 @@
                     vectors = new Vector *[count];
                     j = 0;
-                    first = mol->start;
-                    while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
-                      first = first->next;
-                      Elements[j] = first->type;
-                      vectors[j] = &first->x;
+                    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+                      Elements[j] = (*iter)->type;
+                      vectors[j] = &(*iter)->x;
                       j++;
                     }
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/config.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -1542,5 +1542,4 @@
   int AtomNo = -1;
   int MolNo = 0;
-  atom *Walker = NULL;
   FILE *f = NULL;
 
@@ -1555,26 +1554,24 @@
   fprintf(f, "# Created by MoleCuilder\n");
 
-  for (MoleculeList::const_iterator Runner = MolList->ListOfMolecules.begin(); Runner != MolList->ListOfMolecules.end(); Runner++) {
-    Walker = (*Runner)->start;
+  for (MoleculeList::const_iterator MolRunner = MolList->ListOfMolecules.begin(); MolRunner != MolList->ListOfMolecules.end(); MolRunner++) {
     int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
     AtomNo = 0;
-    while (Walker->next != (*Runner)->end) {
-      Walker = Walker->next;
-      sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
-      elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
+    for (molecule::const_iterator iter = (*MolRunner)->begin(); iter != (*MolRunner)->end(); ++iter) {
+      sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
+      elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
       fprintf(f,
              "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
-             Walker->nr,                /* atom serial number */
+             (*iter)->nr,                /* atom serial number */
              name,         /* atom name */
-             (*Runner)->name,      /* residue name */
+             (*MolRunner)->name,      /* residue name */
              'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
              MolNo,         /* residue sequence number */
-             Walker->node->x[0],                 /* position X in Angstroem */
-             Walker->node->x[1],                 /* position Y in Angstroem */
-             Walker->node->x[2],                 /* position Z in Angstroem */
-             (double)Walker->type->Valence,         /* occupancy */
-             (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
+             (*iter)->node->x[0],                 /* position X in Angstroem */
+             (*iter)->node->x[1],                 /* position Y in Angstroem */
+             (*iter)->node->x[2],                 /* position Z in Angstroem */
+             (double)(*iter)->type->Valence,         /* occupancy */
+             (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
              "0",            /* segment identifier */
-             Walker->type->symbol,    /* element symbol */
+             (*iter)->type->symbol,    /* element symbol */
              "0");           /* charge */
       AtomNo++;
@@ -1596,5 +1593,4 @@
 {
   int AtomNo = -1;
-  atom *Walker = NULL;
   FILE *f = NULL;
 
@@ -1611,24 +1607,22 @@
   fprintf(f, "# Created by MoleCuilder\n");
 
-  Walker = mol->start;
   AtomNo = 0;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
-    elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
+    elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
     fprintf(f,
            "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
-           Walker->nr,                /* atom serial number */
+           (*iter)->nr,                /* atom serial number */
            name,         /* atom name */
            mol->name,      /* residue name */
            'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
            0,         /* residue sequence number */
-           Walker->node->x[0],                 /* position X in Angstroem */
-           Walker->node->x[1],                 /* position Y in Angstroem */
-           Walker->node->x[2],                 /* position Z in Angstroem */
-           (double)Walker->type->Valence,         /* occupancy */
-           (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
+           (*iter)->node->x[0],                 /* position X in Angstroem */
+           (*iter)->node->x[1],                 /* position Y in Angstroem */
+           (*iter)->node->x[2],                 /* position Z in Angstroem */
+           (double)(*iter)->type->Valence,         /* occupancy */
+           (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
            "0",            /* segment identifier */
-           Walker->type->symbol,    /* element symbol */
+           (*iter)->type->symbol,    /* element symbol */
            "0");           /* charge */
     AtomNo++;
@@ -1647,5 +1641,4 @@
 bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
 {
-  atom *Walker = NULL;
   ofstream *output = NULL;
   stringstream * const fname = new stringstream;
@@ -1660,27 +1653,24 @@
 
   // scan maximum number of neighbours
-  Walker = mol->start;
   int MaxNeighbours = 0;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    const int count = Walker->ListOfBonds.size();
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    const int count = (*iter)->ListOfBonds.size();
     if (MaxNeighbours < count)
       MaxNeighbours = count;
   }
-  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
-
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    *output << Walker->nr << "\t";
-    *output << Walker->Name << "\t";
+  *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
+
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    *output << (*iter)->nr << "\t";
+    *output << (*iter)->Name << "\t";
     *output << mol->name << "\t";
     *output << 0 << "\t";
-    *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
-    *output << (double)Walker->type->Valence << "\t";
-    *output << Walker->type->symbol << "\t";
-    for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
-      *output << (*runner)->GetOtherAtom(Walker)->nr << "\t";
-    for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
+    for (int i=0;i<NDIM;i++)
+      *output << (*iter)->node->x[i] << "\t";
+    *output << (double)(*iter)->type->Valence << "\t";
+    *output << (*iter)->type->symbol << "\t";
+    for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
+      *output << (*runner)->GetOtherAtom((*iter))->nr << "\t";
+    for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
       *output << "-\t";
     *output << endl;
@@ -1701,5 +1691,4 @@
 bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
 {
-  atom *Walker = NULL;
   ofstream *output = NULL;
   stringstream * const fname = new stringstream;
@@ -1716,13 +1705,11 @@
   int MaxNeighbours = 0;
   for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
-    Walker = (*MolWalker)->start;
-    while (Walker->next != (*MolWalker)->end) {
-      Walker = Walker->next;
-      const int count = Walker->ListOfBonds.size();
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      const int count = (*iter)->ListOfBonds.size();
       if (MaxNeighbours < count)
         MaxNeighbours = count;
     }
   }
-  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
+  *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
 
   // create global to local id map
@@ -1745,17 +1732,16 @@
     int AtomNo = 0;
     for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
-      Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         *output << AtomNo << "\t";
-        *output << Walker->Name << "\t";
+        *output << (*iter)->Name << "\t";
         *output << (*MolWalker)->name << "\t";
         *output << MolCounter << "\t";
-        *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
-        *output << (double)Walker->type->Valence << "\t";
-        *output << Walker->type->symbol << "\t";
-        for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
-          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ] << "\t";
-        for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
+        for (int i=0;i<NDIM;i++)
+          *output << (*iter)->node->x[i] << "\t";
+        *output << (double)(*iter)->type->Valence << "\t";
+        *output << (*iter)->type->symbol << "\t";
+        for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
+          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom((*iter))->nr ] << "\t";
+        for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
           *output << "-\t";
         *output << endl;
Index: src/helpers.hpp
===================================================================
--- src/helpers.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/helpers.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -100,60 +100,4 @@
 };
 
-/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
- * \param *start begin of chain list
- * \paran *end end of chain list
- * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
- * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
- * \return true - success, false - failure
- */
-template <typename T> bool CreateFatherLookupTable(T *start, T *end, T **&LookupTable, int count = 0)
-{
-  bool status = true;
-  T *Walker = NULL;
-  int AtomNo;
-
-  if (LookupTable != NULL) {
-    Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
-    return false;
-  }
-
-  // count them
-  if (count == 0) {
-    Walker = start;
-    while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
-      Walker = Walker->next;
-      count = (count < Walker->GetTrueFather()->nr) ? Walker->GetTrueFather()->nr : count;
-    }
-  }
-  if (count <= 0) {
-    Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
-    return false;
-  }
-
-  // allocate and fill
-  LookupTable = Calloc<T*>(count, "CreateFatherLookupTable - **LookupTable");
-  if (LookupTable == NULL) {
-    eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
-    performCriticalExit();
-    status = false;
-  } else {
-    Walker = start;
-    while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
-      Walker = Walker->next;
-      AtomNo = Walker->GetTrueFather()->nr;
-      if ((AtomNo >= 0) && (AtomNo < count)) {
-        //*out << "Setting LookupTable[" << AtomNo << "] to " << *Walker << endl;
-        LookupTable[AtomNo] = Walker;
-      } else {
-        Log() << Verbose(0) << "Walker " << *Walker << " exceeded range of nuclear ids [0, " << count << ")." << endl;
-        status = false;
-        break;
-      }
-    }
-  }
-
-  return status;
-};
-
 /** Frees a two-dimensional array.
  * \param *ptr pointer to array
Index: src/lists.cpp
===================================================================
--- src/lists.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ 	(revision )
@@ -1,21 +1,0 @@
-/*
- * lists.cpp
- *
- *  Created on: Feb 24, 2010
- *      Author: crueger
- */
-
-#include "lists.hpp"
-#include "World.hpp"
-#include "atom.hpp"
-
-// specialized templates for lists
-// can be directly compiled since no further templatization is needed
-template <> void removewithoutcheck<atom>(atom *walker)
-{
-  if (walker != 0) {
-    unlink(walker);
-    walker->getWorld()->destroyAtom(walker);
-    walker = 0;
-  }
-};
Index: src/lists.hpp
===================================================================
--- src/lists.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/lists.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -8,6 +8,4 @@
 #ifndef LISTS_HPP_
 #define LISTS_HPP_
-
-class atom;
 
 /******************************** Some templates for list management ***********************************/
@@ -134,27 +132,3 @@
 };
 
-/** Returns the first marker in a chain list.
- * \param *me one arbitrary item in chain list
- * \return poiner to first marker
- */
-template <typename X> X *GetFirst(X *me)
-{
-  X *Binder = me;
-  while(Binder->previous != 0)
-    Binder = Binder->previous;
-  return Binder;
-};
-
-/** Returns the last marker in a chain list.
- * \param *me one arbitrary item in chain list
- * \return poiner to last marker
- */
-template <typename X> X *GetLast(X *me)
-{
-  X *Binder = me;
-  while(Binder->next != 0)
-    Binder = Binder->next;
-  return Binder;
-};
-
 #endif /* LISTS_HPP_ */
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -37,11 +37,6 @@
   formula(this,boost::bind(&molecule::calcFormula,this)),
   last_atom(0),
-  InternalPointer(start)
-{
-  // init atom chain list
-  //start->father = NULL;
-  //end->father = NULL;
-  //link(start,end);
-
+  InternalPointer(begin())
+{
   // init bond chain list
   link(first,last);
@@ -67,6 +62,4 @@
   delete(first);
   delete(last);
-  end->getWorld()->destroyAtom(end);
-  start->getWorld()->destroyAtom(start);
 };
 
@@ -103,6 +96,6 @@
   for (int j = 0; j<MAX_ELEMENTS;j++)
     Counts[j] = 0;
-  for(atom *Walker = start; Walker != end; Walker = Walker->next) {
-    Counts[Walker->type->Z]++;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    Counts[(*iter)->type->Z]++;
   }
   for(element* Elemental = elemente->end; Elemental != elemente->start; Elemental = Elemental->previous) {
@@ -124,10 +117,54 @@
 }
 
-molecule::iterator molecule::ende(){
+molecule::iterator molecule::end(){
   return molecule::iterator(atoms.end(),this);
 }
 
-molecule::const_iterator molecule::ende() const{
+molecule::const_iterator molecule::end() const{
   return atoms.end();
+}
+
+bool molecule::empty() const
+{
+  return (begin() == end());
+}
+
+size_t molecule::size() const
+{
+  size_t counter = 0;
+  for (molecule::const_iterator iter = begin(); iter != end (); ++iter)
+    counter++;
+  return counter;
+}
+
+molecule::const_iterator molecule::erase( const_iterator loc )
+{
+  molecule::const_iterator iter = loc;
+  iter--;
+  atoms.erase( loc );
+  return iter;
+}
+
+molecule::const_iterator molecule::erase( atom *& key )
+{
+  molecule::const_iterator iter = find(key);
+  molecule::const_iterator runner = iter;
+  if (runner != begin()) {
+    runner--;
+    if (iter != end())
+      atoms.erase( key );
+    return runner;
+  } else
+    return end();
+}
+
+molecule::const_iterator molecule::find ( atom *& key ) const
+{
+  return atoms.find( key );
+}
+
+pair<molecule::iterator,bool> molecule::insert ( atom * const key )
+{
+  return atoms.insert(key);
 }
 
@@ -139,9 +176,7 @@
 bool molecule::AddAtom(atom *pointer)
 {
-  bool retval = false;
   OBSERVE;
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    pointer->nr = last_atom++;  // increase number within molecule
     AtomCount++;
     if (pointer->type != NULL) {
@@ -157,7 +192,7 @@
       }
     }
-    retval = add(pointer, end);
-  }
-  return retval;
+    insert(pointer);
+  }
+  return true;
 };
 
@@ -175,6 +210,5 @@
     walker->Name = Malloc<char>(strlen(pointer->Name) + 1, "atom::atom: *Name");
     strcpy (walker->Name, pointer->Name);
-    walker->nr = last_atom++;  // increase number within molecule
-    add(walker, end);
+    insert(walker);
     if ((pointer->type != NULL) && (pointer->type->Z != 1))
       NoNonHydrogen++;
@@ -701,5 +735,6 @@
     ElementCount--;
   RemoveBonds(pointer);
-  return remove(pointer, start, end);
+  erase(pointer);
+  return true;
 };
 
@@ -718,5 +753,5 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
-  unlink(pointer);
+  erase(pointer);
   return true;
 };
@@ -727,5 +762,7 @@
 bool molecule::CleanupMolecule()
 {
-  return (cleanup(first,last) && cleanup(start,end));
+  for (molecule::iterator iter = begin(); !empty(); iter = begin())
+      erase(iter);
+  return (cleanup(first,last));
 };
 
@@ -734,9 +771,13 @@
  * \return pointer to atom or NULL
  */
-atom * molecule::FindAtom(int Nr)  const{
-  atom * walker = find(&Nr, start,end);
-  if (walker != NULL) {
+atom * molecule::FindAtom(int Nr)  const
+{
+  molecule::const_iterator iter = begin();
+  for (; iter != end(); ++iter)
+    if ((*iter)->nr == Nr)
+      break;
+  if (iter != end()) {
     //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
-    return walker;
+    return (*iter);
   } else {
     Log() << Verbose(0) << "Atom not found in list." << endl;
@@ -896,10 +937,5 @@
 void molecule::CountAtoms()
 {
-  int i = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    i++;
-  }
+  int i = size();
   if ((AtomCount == 0) || (i != AtomCount)) {
     Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
@@ -910,14 +946,12 @@
       i=0;
       NoNonHydrogen = 0;
-      Walker = start;
-      while (Walker->next != end) {
-        Walker = Walker->next;
-        Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-        if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
+      for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+        (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
+        if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
           NoNonHydrogen++;
-        Free(&Walker->Name);
-        Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
-        sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
-        Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
+        Free(&(*iter)->Name);
+        (*iter)->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
+        sprintf((*iter)->Name, "%2s%02d", (*iter)->type->symbol, (*iter)->nr+1);
+        Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->Name << "." << endl;
         i++;
       }
@@ -1092,5 +1126,4 @@
 int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
 {
-  atom *Walker = NULL, *OtherWalker = NULL;
   Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
   int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
@@ -1103,21 +1136,17 @@
   } else {
     Log() << Verbose(4) << "Map is ";
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      if (Walker->father == NULL) {
-        AtomicMap[Walker->nr] = -2;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+      if ((*iter)->father == NULL) {
+        AtomicMap[(*iter)->nr] = -2;
       } else {
-        OtherWalker = OtherMolecule->start;
-        while (OtherWalker->next != OtherMolecule->end) {
-          OtherWalker = OtherWalker->next;
+        for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
       //for (int i=0;i<AtomCount;i++) { // search atom
         //for (int j=0;j<OtherMolecule->AtomCount;j++) {
-          //Log() << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
-          if (Walker->father == OtherWalker)
-            AtomicMap[Walker->nr] = OtherWalker->nr;
+          //Log() << Verbose(4) << "Comparing father " << (*iter)->father << " with the other one " << (*runner)->father << "." << endl;
+          if ((*iter)->father == (*runner))
+            AtomicMap[(*iter)->nr] = (*runner)->nr;
         }
       }
-      Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";
+      Log() << Verbose(0) << AtomicMap[(*iter)->nr] << "\t";
     }
     Log() << Verbose(0) << endl;
@@ -1153,8 +1182,6 @@
 void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    array[(Walker->*index)] = Walker;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    array[((*iter)->*index)] = (*iter);
   }
 };
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -137,4 +137,10 @@
   iterator end();
   const_iterator end() const;
+  bool empty() const;
+  size_t size() const;
+  const_iterator erase( const_iterator loc );
+  const_iterator erase( atom *& key );
+  const_iterator find (  atom *& key ) const;
+  pair<iterator,bool> insert ( atom * const key );
 
 
@@ -143,10 +149,7 @@
   Vector *GetCenter() const ;
   TesselPoint *GetPoint() const ;
-  TesselPoint *GetTerminalPoint() const ;
   int GetMaxId() const;
   void GoToNext() const ;
-  void GoToPrevious() const ;
   void GoToFirst() const ;
-  void GoToLast() const ;
   bool IsEmpty() const ;
   bool IsEnd() const ;
@@ -310,4 +313,5 @@
   bool StoreForcesFile(MoleculeListClass *BondFragments, char *path, int *SortIndex);
   bool CreateMappingLabelsToConfigSequence(int *&SortIndex);
+  bool CreateFatherLookupTable(atom **&LookupTable, int count = 0);
   void BreadthFirstSearchAdd(molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
   /// -# BOSSANOVA
@@ -338,5 +342,5 @@
   private:
   int last_atom;      //!< number given to last atom
-  mutable atom *InternalPointer;  //!< internal pointer for PointCloud
+  mutable iterator InternalPointer;  //!< internal pointer for PointCloud
 };
 
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule_dynamics.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -27,12 +27,10 @@
   gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
   gsl_vector *x = gsl_vector_alloc(NDIM);
-  atom * Runner = mol->start;
   atom *Sprinter = NULL;
   Vector trajectory1, trajectory2, normal, TestVector;
   double Norm1, Norm2, tmp, result = 0.;
 
-  while (Runner->next != mol->end) {
-    Runner = Runner->next;
-    if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
       break;
     // determine normalized trajectories direction vector (n1, n2)
@@ -42,21 +40,21 @@
     trajectory1.Normalize();
     Norm1 = trajectory1.Norm();
-    Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
+    Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
     trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
-    trajectory2.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));
+    trajectory2.SubtractVector(&(*iter)->Trajectory.R.at(Params.startstep));
     trajectory2.Normalize();
     Norm2 = trajectory1.Norm();
     // check whether either is zero()
     if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
-      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
+      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&(*iter)->Trajectory.R.at(Params.startstep));
     } else if (Norm1 < MYEPSILON) {
       Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
       trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy first offset
-      trajectory1.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));  // subtract second offset
+      trajectory1.SubtractVector(&(*iter)->Trajectory.R.at(Params.startstep));  // subtract second offset
       trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
       trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
       tmp = trajectory1.Norm();  // remaining norm is distance
     } else if (Norm2 < MYEPSILON) {
-      Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
+      Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
       trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy second offset
       trajectory2.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));  // subtract first offset
@@ -65,12 +63,12 @@
       tmp = trajectory2.Norm();  // remaining norm is distance
     } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
-  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
+  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear dependent: ";
   //        Log() << Verbose(0) << trajectory1;
   //        Log() << Verbose(0) << " and ";
   //        Log() << Verbose(0) << trajectory2;
-      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
+      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&(*iter)->Trajectory.R.at(Params.startstep));
   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
     } else { // determine distance by finding minimum distance
-  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
+  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
   //        Log() << Verbose(0) << endl;
   //        Log() << Verbose(0) << "First Trajectory: ";
@@ -88,5 +86,5 @@
         gsl_matrix_set(A, 1, i, trajectory2.x[i]);
         gsl_matrix_set(A, 2, i, normal.x[i]);
-        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep).x[i] - Runner->Trajectory.R.at(Params.startstep).x[i]));
+        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep).x[i] - (*iter)->Trajectory.R.at(Params.startstep).x[i]));
       }
       // solve the linear system by Householder transformations
@@ -96,5 +94,5 @@
   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
       // test whether we really have the intersection (by checking on c_1 and c_2)
-      TestVector.CopyVector(&Runner->Trajectory.R.at(Params.startstep));
+      TestVector.CopyVector(&(*iter)->Trajectory.R.at(Params.startstep));
       trajectory2.Scale(gsl_vector_get(x,1));
       TestVector.AddVector(&trajectory2);
@@ -131,10 +129,8 @@
 {
   double result = 0.;
-  atom * Runner = mol->start;
-  while (Runner->next != mol->end) {
-    Runner = Runner->next;
-    if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[(*iter)->nr]) && (Walker->nr < (*iter)->nr)) {
   //    atom *Sprinter = PermutationMap[Walker->nr];
-  //        Log() << Verbose(0) << *Walker << " and " << *Runner << " are heading to the same target at ";
+  //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
   //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
   //        Log() << Verbose(0) << ", penalting." << endl;
@@ -167,10 +163,8 @@
   // go through every atom
   atom *Runner = NULL;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     // first term: distance to target
-    Runner = Params.PermutationMap[Walker->nr];   // find target point
-    tmp = (Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)));
+    Runner = Params.PermutationMap[(*iter)->nr];   // find target point
+    tmp = ((*iter)->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)));
     tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     result += Params.PenaltyConstants[0] * tmp;
@@ -178,8 +172,8 @@
 
     // second term: sum of distances to other trajectories
-    result += SumDistanceOfTrajectories(Walker, this, Params);
+    result += SumDistanceOfTrajectories((*iter), this, Params);
 
     // third term: penalty for equal targets
-    result += PenalizeEqualTargets(Walker, this, Params);
+    result += PenalizeEqualTargets((*iter), this, Params);
   }
 
@@ -224,12 +218,7 @@
   }
 
-  atom *Runner = NULL;
-  atom *Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    Runner = mol->start;
-    while(Runner->next != mol->end) {
-      Runner = Runner->next;
-      Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)), Runner) );
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
+      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).Distance(&(*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
     }
   }
@@ -243,12 +232,10 @@
 void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
 {
-  atom *Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    Params.StepList[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
-    Params.PermutationMap[Walker->nr] = Params.DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
-    Params.DoubleList[Params.DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
-    Params.DistanceIterators[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // and remember which one we picked
-    Log() << Verbose(2) << *Walker << " starts with distance " << Params.DistanceList[Walker->nr]->begin()->first << "." << endl;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    Params.StepList[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // stores the step to the next iterator that could be a possible next target
+    Params.PermutationMap[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin()->second;   // always pick target with the smallest distance
+    Params.DoubleList[Params.DistanceList[(*iter)->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
+    Params.DistanceIterators[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // and remember which one we picked
+    Log() << Verbose(2) << *(*iter) << " starts with distance " << Params.DistanceList[(*iter)->nr]->begin()->first << "." << endl;
   }
 };
@@ -291,17 +278,21 @@
 void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
 {
-  atom *Walker = mol->start;
+  molecule::const_iterator iter = mol->begin();
   DistanceMap::iterator NewBase;
   double Potential = fabs(mol->ConstrainedPotential(Params));
 
+  if (mol->empty()) {
+    eLog() << Verbose(1) << "Molecule is empty." << endl;
+    return;
+  }
   while ((Potential) > Params.PenaltyConstants[2]) {
     PrintPermutationMap(mol->AtomCount, Params);
-    Walker = Walker->next;
-    if (Walker == mol->end) // round-robin at the end
-      Walker = mol->start->next;
-    if (Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
+    iter++;
+    if (iter == mol->end()) // round-robin at the end
+      iter = mol->begin();
+    if (Params.DoubleList[Params.DistanceIterators[(*iter)->nr]->second->nr] <= 1)  // no need to make those injective that aren't
       continue;
     // now, try finding a new one
-    Potential = TryNextNearestNeighbourForInjectivePermutation(mol, Walker, Potential, Params);
+    Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params);
   }
   for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1
@@ -350,5 +341,5 @@
   Params.StepList = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*StepList");
   int round;
-  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
+  atom *Sprinter = NULL;
   DistanceMap::iterator Rider, Strider;
 
@@ -377,37 +368,36 @@
     Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl;
     OlderPotential = OldPotential;
+    molecule::const_iterator iter;
     do {
-      Walker = start;
-      while (Walker->next != end) { // pick one
-        Walker = Walker->next;
+      iter = begin();
+      for (; iter != end(); ++iter) {
         PrintPermutationMap(AtomCount, Params);
-        Sprinter = Params.DistanceIterators[Walker->nr]->second;   // store initial partner
-        Strider = Params.DistanceIterators[Walker->nr];  //remember old iterator
-        Params.DistanceIterators[Walker->nr] = Params.StepList[Walker->nr];
-        if (Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
-          Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->begin();
+        Sprinter = Params.DistanceIterators[(*iter)->nr]->second;   // store initial partner
+        Strider = Params.DistanceIterators[(*iter)->nr];  //remember old iterator
+        Params.DistanceIterators[(*iter)->nr] = Params.StepList[(*iter)->nr];
+        if (Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->end()) {// stop, before we run through the list and still on
+          Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->begin();
           break;
         }
-        //Log() << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
+        //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->nr]->second << "." << endl;
         // find source of the new target
-        Runner = start->next;
-        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
-          if (Params.PermutationMap[Runner->nr] == Params.DistanceIterators[Walker->nr]->second) {
-            //Log() << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
+        molecule::const_iterator runner = begin();
+        for (; runner != end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
+          if (Params.PermutationMap[(*runner)->nr] == Params.DistanceIterators[(*iter)->nr]->second) {
+            //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->nr] << "." << endl;
             break;
           }
-          Runner = Runner->next;
         }
-        if (Runner != end) { // we found the other source
+        if (runner != end()) { // we found the other source
           // then look in its distance list for Sprinter
-          Rider = Params.DistanceList[Runner->nr]->begin();
-          for (; Rider != Params.DistanceList[Runner->nr]->end(); Rider++)
+          Rider = Params.DistanceList[(*runner)->nr]->begin();
+          for (; Rider != Params.DistanceList[(*runner)->nr]->end(); Rider++)
             if (Rider->second == Sprinter)
               break;
-          if (Rider != Params.DistanceList[Runner->nr]->end()) { // if we have found one
-            //Log() << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
+          if (Rider != Params.DistanceList[(*runner)->nr]->end()) { // if we have found one
+            //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->nr] << "/" << *Rider->second << "." << endl;
             // exchange both
-            Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
-            Params.PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
+            Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // put next farther distance into PermutationMap
+            Params.PermutationMap[(*runner)->nr] = Sprinter;  // and hand the old target to its respective owner
             PrintPermutationMap(AtomCount, Params);
             // calculate the new potential
@@ -416,13 +406,13 @@
             if (Potential > OldPotential) { // we made everything worse! Undo ...
               //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
-              //Log() << Verbose(3) << "Setting " << *Runner << "'s source to " << *Params.DistanceIterators[Runner->nr]->second << "." << endl;
+              //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->nr]->second << "." << endl;
               // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
-              Params.PermutationMap[Runner->nr] = Params.DistanceIterators[Runner->nr]->second;
+              Params.PermutationMap[(*runner)->nr] = Params.DistanceIterators[(*runner)->nr]->second;
               // Undo for Walker
-              Params.DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
-              //Log() << Verbose(3) << "Setting " << *Walker << "'s source to " << *Params.DistanceIterators[Walker->nr]->second << "." << endl;
-              Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
+              Params.DistanceIterators[(*iter)->nr] = Strider;  // take next farther distance target
+              //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->nr]->second << "." << endl;
+              Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second;
             } else {
-              Params.DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
+              Params.DistanceIterators[(*runner)->nr] = Rider;  // if successful also move the pointer in the iterator list
               Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
               OldPotential = Potential;
@@ -434,13 +424,13 @@
             //Log() << Verbose(0) << endl;
           } else {
-            eLog() << Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
+            eLog() << Verbose(1) << *(*runner) << " was not the owner of " << *Sprinter << "!" << endl;
             exit(255);
           }
         } else {
-          Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // new target has no source!
+          Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // new target has no source!
         }
-        Params.StepList[Walker->nr]++; // take next farther distance target
+        Params.StepList[(*iter)->nr]++; // take next farther distance target
       }
-    } while (Walker->next != end);
+    } while (++iter != end());
   } while ((OlderPotential - OldPotential) > 1e-3);
   Log() << Verbose(1) << "done." << endl;
@@ -489,5 +479,5 @@
   // Get the Permutation Map by MinimiseConstrainedPotential
   atom **PermutationMap = NULL;
-  atom *Walker = NULL, *Sprinter = NULL;
+  atom *Sprinter = NULL;
   if (!MapByIdentity)
     MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
@@ -508,17 +498,15 @@
     mol = World::get()->createMolecule();
     MoleculePerStep->insert(mol);
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
       // add to molecule list
-      Sprinter = mol->AddCopyAtom(Walker);
+      Sprinter = mol->AddCopyAtom((*iter));
       for (int n=NDIM;n--;) {
-        Sprinter->x.x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+        Sprinter->x.x[n] = (*iter)->Trajectory.R.at(startstep).x[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep).x[n] - (*iter)->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
         // add to Trajectories
         //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
         if (step < MaxSteps) {
-          Walker->Trajectory.R.at(step).x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
-          Walker->Trajectory.U.at(step).x[n] = 0.;
-          Walker->Trajectory.F.at(step).x[n] = 0.;
+          (*iter)->Trajectory.R.at(step).x[n] = (*iter)->Trajectory.R.at(startstep).x[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep).x[n] - (*iter)->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+          (*iter)->Trajectory.U.at(step).x[n] = 0.;
+          (*iter)->Trajectory.F.at(step).x[n] = 0.;
         }
       }
Index: src/molecule_fragmentation.cpp
===================================================================
--- src/molecule_fragmentation.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule_fragmentation.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -38,8 +38,6 @@
   int FragmentCount;
   // get maximum bond degree
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    c = (Walker->ListOfBonds.size() > c) ? Walker->ListOfBonds.size() : c;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c;
   }
   FragmentCount = NoNonHydrogen*(1 << (c*order));
@@ -358,5 +356,5 @@
 map<double, pair<int,int> >  * ReMapAdaptiveCriteriaListToValue(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
 {
-  atom *Walker = mol->start;
+  atom *Walker = NULL;
   map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
   Log() << Verbose(1) << "Root candidate list is: " << endl;
@@ -390,5 +388,5 @@
 bool MarkUpdateCandidates(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
 {
-  atom *Walker = mol->start;
+  atom *Walker = NULL;
   int No = -1;
   bool status = false;
@@ -434,5 +432,4 @@
 bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
 {
-  atom *Walker = start;
   bool status = false;
 
@@ -456,11 +453,10 @@
     if (AdaptiveCriteriaList->empty()) {
       eLog() << Verbose(2) << "Unable to parse file, incrementing all." << endl;
-      while (Walker->next != end) {
-        Walker = Walker->next;
+      for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     #ifdef ADDHYDROGEN
-        if (Walker->type->Z != 1) // skip hydrogen
+        if ((*iter)->type->Z != 1) // skip hydrogen
     #endif
         {
-          AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
+          AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
           status = true;
         }
@@ -477,12 +473,11 @@
     Free(&FinalRootCandidates);
   } else { // global increase of Bond Order
-    while (Walker->next != end) {
-      Walker = Walker->next;
+    for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
   #ifdef ADDHYDROGEN
-      if (Walker->type->Z != 1) // skip hydrogen
+      if ((*iter)->type->Z != 1) // skip hydrogen
   #endif
       {
-        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
-        if ((Order != 0) && (Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
+        AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
+        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->nr]))
           status = true;
       }
@@ -523,4 +518,57 @@
 
   return true;
+};
+
+
+
+/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
+ * \param *start begin of list (STL iterator, i.e. first item)
+ * \paran *end end of list (STL iterator, i.e. one past last item)
+ * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
+ * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
+ * \return true - success, false - failure
+ */
+bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
+{
+  bool status = true;
+  int AtomNo;
+
+  if (LookupTable != NULL) {
+    Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
+    return false;
+  }
+
+  // count them
+  if (count == 0) {
+    for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
+      count = (count < (*iter)->GetTrueFather()->nr) ? (*iter)->GetTrueFather()->nr : count;
+    }
+  }
+  if (count <= 0) {
+    Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
+    return false;
+  }
+
+  // allocate and fill
+  LookupTable = Calloc<atom *>(count, "CreateFatherLookupTable - **LookupTable");
+  if (LookupTable == NULL) {
+    eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
+    performCriticalExit();
+    status = false;
+  } else {
+    for (molecule::iterator iter = begin(); iter != end(); ++iter) {
+      AtomNo = (*iter)->GetTrueFather()->nr;
+      if ((AtomNo >= 0) && (AtomNo < count)) {
+        //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl;
+        LookupTable[AtomNo] = (*iter);
+      } else {
+        Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl;
+        status = false;
+        break;
+      }
+    }
+  }
+
+  return status;
 };
 
@@ -575,5 +623,5 @@
 
   // create lookup table for Atom::nr
-  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(start, end, ListOfAtoms, AtomCount);
+  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, AtomCount);
 
   // === compare it with adjacency file ===
@@ -871,9 +919,10 @@
   Leaf->CountAtoms();
 
-  atom *Runner = Leaf->start;
-  while (Runner->next != Leaf->end) {
-    Runner = Runner->next;
+#ifdef ADDHYDROGEN
+  molecule::const_iterator runner;
+#endif
+  for (molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ) {
     LonelyFlag = true;
-    FatherOfRunner = Runner->father;
+    FatherOfRunner = (*iter)->father;
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
@@ -886,7 +935,7 @@
 //            Log() << Verbose(3) << "Adding Bond: ";
 //            Log() << Verbose(0) <<
-            Leaf->AddBond(Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree);
+            Leaf->AddBond((*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree);
 //            Log() << Verbose(0) << "." << endl;
-            //NumBonds[Runner->nr]++;
+            //NumBonds[(*iter)->nr]++;
           } else {
 //            Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
@@ -896,20 +945,20 @@
 //          Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl;
 #ifdef ADDHYDROGEN
-          //Log() << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), Runner, FatherOfRunner, OtherFather, IsAngstroem))
+          //Log() << Verbose(3) << "Adding Hydrogen to " << (*iter)->Name << " and a bond in between." << endl;
+          if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
             exit(1);
 #endif
-          //NumBonds[Runner->nr] += Binder->BondDegree;
+          //NumBonds[(*iter)->nr] += Binder->BondDegree;
         }
       }
     } else {
-      eLog() << Verbose(1) << "Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
+      eLog() << Verbose(1) << "Son " << (*iter)->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
     }
     if ((LonelyFlag) && (Leaf->AtomCount > 1)) {
-      Log() << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl;
+      Log() << Verbose(0) << *(*iter) << "has got bonds only to hydrogens!" << endl;
     }
 #ifdef ADDHYDROGEN
-    while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
-      Runner = Runner->next;
+    while ((iter != Leaf->end()) && ((*iter)->type->Z == 1)) // skip added hydrogen
+      iter++;
 #endif
   }
@@ -1559,9 +1608,7 @@
 
   // Construct the complete KeySet which we need for topmost level only (but for all Roots)
-  atom *Walker = start;
   KeySet CompleteMolecule;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    CompleteMolecule.insert(Walker->GetTrueFather()->nr);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    CompleteMolecule.insert((*iter)->GetTrueFather()->nr);
   }
 
@@ -1574,5 +1621,5 @@
     RootKeyNr = RootStack.front();
     RootStack.pop_front();
-    Walker = FindAtom(RootKeyNr);
+    atom *Walker = FindAtom(RootKeyNr);
     // check cyclic lengths
     //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule_geometry.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -66,16 +66,15 @@
 
 //  Log() << Verbose(3) << "Begin of CenterEdge." << endl;
-  atom *ptr = start->next;  // start at first in list
-  if (ptr != end) {   //list not empty?
+  molecule::const_iterator iter = begin();  // start at first in list
+  if (iter != end()) { //list not empty?
     for (int i=NDIM;i--;) {
-      max->x[i] = ptr->x.x[i];
-      min->x[i] = ptr->x.x[i];
-    }
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
-      //ptr->Output(1,1,out);
+      max->x[i] = (*iter)->x.x[i];
+      min->x[i] = (*iter)->x.x[i];
+    }
+    for (; iter != end(); ++iter) {// continue with second if present
+      //(*iter)->Output(1,1,out);
       for (int i=NDIM;i--;) {
-        max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
-        min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i];
+        max->x[i] = (max->x[i] < (*iter)->x.x[i]) ? (*iter)->x.x[i] : max->x[i];
+        min->x[i] = (min->x[i] > (*iter)->x.x[i]) ? (*iter)->x.x[i] : min->x[i];
       }
     }
@@ -101,13 +100,12 @@
 {
   int Num = 0;
-  atom *ptr = start;  // start at first in list
+  molecule::const_iterator iter = begin();  // start at first in list
 
   Center.Zero();
 
-  if (ptr->next != end) {   //list not empty?
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
+  if (iter != end()) {   //list not empty?
+    for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      Center.AddVector(&ptr->x);
+      Center.AddVector(&(*iter)->x);
     }
     Center.Scale(-1./Num); // divide through total number (and sign for direction)
@@ -122,5 +120,5 @@
 Vector * molecule::DetermineCenterOfAll() const
 {
-  atom *ptr = start->next;  // start at first in list
+  molecule::const_iterator iter = begin();  // start at first in list
   Vector *a = new Vector();
   Vector tmp;
@@ -129,9 +127,8 @@
   a->Zero();
 
-  if (ptr != end) {   //list not empty?
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
+  if (iter != end()) {   //list not empty?
+    for (; iter != end(); ++iter) {  // continue with second if present
       Num += 1.;
-      tmp.CopyVector(&ptr->x);
+      tmp.CopyVector(&(*iter)->x);
       a->AddVector(&tmp);
     }
@@ -147,5 +144,5 @@
 Vector * molecule::DetermineCenterOfGravity()
 {
-  atom *ptr = start->next;  // start at first in list
+  molecule::const_iterator iter = begin();  // start at first in list
   Vector *a = new Vector();
   Vector tmp;
@@ -154,10 +151,9 @@
   a->Zero();
 
-  if (ptr != end) {   //list not empty?
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
-      Num += ptr->type->mass;
-      tmp.CopyVector(&ptr->x);
-      tmp.Scale(ptr->type->mass);  // scale by mass
+  if (iter != end()) {   //list not empty?
+    for (; iter != end(); ++iter) {  // continue with second if present
+      Num += (*iter)->type->mass;
+      tmp.CopyVector(&(*iter)->x);
+      tmp.Scale((*iter)->type->mass);  // scale by mass
       a->AddVector(&tmp);
     }
@@ -195,11 +191,8 @@
 void molecule::Scale(const double ** const factor)
 {
-  atom *ptr = start;
-
-  while (ptr->next != end) {
-    ptr = ptr->next;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Scale(factor);
-    ptr->x.Scale(factor);
+      (*iter)->Trajectory.R.at(j).Scale(factor);
+    (*iter)->x.Scale(factor);
   }
 };
@@ -210,11 +203,8 @@
 void molecule::Translate(const Vector *trans)
 {
-  atom *ptr = start;
-
-  while (ptr->next != end) {
-    ptr = ptr->next;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Translate(trans);
-    ptr->x.Translate(trans);
+      (*iter)->Trajectory.R.at(j).Translate(trans);
+    (*iter)->x.Translate(trans);
   }
 };
@@ -251,5 +241,4 @@
 void molecule::DeterminePeriodicCenter(Vector &center)
 {
-  atom *Walker = start;
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   double *inversematrix = InverseMatrix(cell_size);
@@ -261,19 +250,18 @@
     Center.Zero();
     flag = true;
-    while (Walker->next != end) {
-      Walker = Walker->next;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 #ifdef ADDHYDROGEN
-      if (Walker->type->Z != 1) {
+      if ((*iter)->type->Z != 1) {
 #endif
-        Testvector.CopyVector(&Walker->x);
+        Testvector.CopyVector(&(*iter)->x);
         Testvector.MatrixMultiplication(inversematrix);
         Translationvector.Zero();
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
-         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
+        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
+         if ((*iter)->nr < (*Runner)->GetOtherAtom((*iter))->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
+              tmp = (*iter)->x.x[j] - (*Runner)->GetOtherAtom((*iter))->x.x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
-                Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl;
+                Log() << Verbose(0) << "Hit: atom " << (*iter)->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl;
                 if (tmp > 0)
                   Translationvector.x[j] -= 1.;
@@ -291,7 +279,7 @@
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
-          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
-            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
+        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
+            Testvector.CopyVector(&(*Runner)->GetOtherAtom((*iter))->x);
             Testvector.MatrixMultiplication(inversematrix);
             Testvector.AddVector(&Translationvector);
@@ -320,5 +308,4 @@
 void molecule::PrincipalAxisSystem(bool DoRotate)
 {
-  atom *ptr = start;  // start at first in list
   double InertiaTensor[NDIM*NDIM];
   Vector *CenterOfGravity = DetermineCenterOfGravity();
@@ -331,18 +318,17 @@
 
   // sum up inertia tensor
-  while (ptr->next != end) {
-    ptr = ptr->next;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     Vector x;
-    x.CopyVector(&ptr->x);
+    x.CopyVector(&(*iter)->x);
     //x.SubtractVector(CenterOfGravity);
-    InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-    InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-    InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-    InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-    InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-    InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-    InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-    InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-    InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+    InertiaTensor[0] += (*iter)->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
+    InertiaTensor[1] += (*iter)->type->mass*(-x.x[0]*x.x[1]);
+    InertiaTensor[2] += (*iter)->type->mass*(-x.x[0]*x.x[2]);
+    InertiaTensor[3] += (*iter)->type->mass*(-x.x[1]*x.x[0]);
+    InertiaTensor[4] += (*iter)->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
+    InertiaTensor[5] += (*iter)->type->mass*(-x.x[1]*x.x[2]);
+    InertiaTensor[6] += (*iter)->type->mass*(-x.x[2]*x.x[0]);
+    InertiaTensor[7] += (*iter)->type->mass*(-x.x[2]*x.x[1]);
+    InertiaTensor[8] += (*iter)->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
   }
   // print InertiaTensor for debugging
@@ -382,19 +368,17 @@
 
     // sum up inertia tensor
-    ptr = start;
-    while (ptr->next != end) {
-      ptr = ptr->next;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
       Vector x;
-      x.CopyVector(&ptr->x);
+      x.CopyVector(&(*iter)->x);
       //x.SubtractVector(CenterOfGravity);
-      InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-      InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-      InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-      InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-      InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-      InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-      InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-      InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-      InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+      InertiaTensor[0] += (*iter)->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
+      InertiaTensor[1] += (*iter)->type->mass*(-x.x[0]*x.x[1]);
+      InertiaTensor[2] += (*iter)->type->mass*(-x.x[0]*x.x[2]);
+      InertiaTensor[3] += (*iter)->type->mass*(-x.x[1]*x.x[0]);
+      InertiaTensor[4] += (*iter)->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
+      InertiaTensor[5] += (*iter)->type->mass*(-x.x[1]*x.x[2]);
+      InertiaTensor[6] += (*iter)->type->mass*(-x.x[2]*x.x[0]);
+      InertiaTensor[7] += (*iter)->type->mass*(-x.x[2]*x.x[1]);
+      InertiaTensor[8] += (*iter)->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
     }
     // print InertiaTensor for debugging
@@ -420,5 +404,4 @@
 void molecule::Align(Vector *n)
 {
-  atom *ptr = start;
   double alpha, tmp;
   Vector z_axis;
@@ -431,13 +414,12 @@
   alpha = atan(-n->x[0]/n->x[2]);
   Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ";
-  while (ptr->next != end) {
-    ptr = ptr->next;
-    tmp = ptr->x.x[0];
-    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
-    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    tmp = (*iter)->x.x[0];
+    (*iter)->x.x[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x.x[2];
+    (*iter)->x.x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x.x[2];
     for (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j).x[0];
-      ptr->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
-      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
+      tmp = (*iter)->Trajectory.R.at(j).x[0];
+      (*iter)->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j).x[2];
+      (*iter)->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j).x[2];
     }
   }
@@ -451,16 +433,14 @@
 
   // rotate on z-y plane
-  ptr = start;
   alpha = atan(-n->x[1]/n->x[2]);
   Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ";
-  while (ptr->next != end) {
-    ptr = ptr->next;
-    tmp = ptr->x.x[1];
-    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
-    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    tmp = (*iter)->x.x[1];
+    (*iter)->x.x[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x.x[2];
+    (*iter)->x.x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x.x[2];
     for (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j).x[1];
-      ptr->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
-      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
+      tmp = (*iter)->Trajectory.R.at(j).x[1];
+      (*iter)->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j).x[2];
+      (*iter)->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j).x[2];
     }
   }
@@ -487,5 +467,4 @@
   Vector a,b,c,d;
   struct lsq_params *par = (struct lsq_params *)params;
-  atom *ptr = par->mol->start;
 
   // initialize vectors
@@ -497,8 +476,7 @@
   b.x[2] = gsl_vector_get(x,5);
   // go through all atoms
-  while (ptr != par->mol->end) {
-    ptr = ptr->next;
-    if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
-      c.CopyVector(&ptr->x);  // copy vector to temporary one
+  for (molecule::const_iterator iter = par->mol->begin(); iter != par->mol->end(); ++iter) {
+    if ((*iter)->type == ((struct lsq_params *)params)->type) { // for specific type
+      c.CopyVector(&(*iter)->x);  // copy vector to temporary one
       c.SubtractVector(&a);   // subtract offset vector
       t = c.ScalarProduct(&b);           // get direction parameter
Index: src/molecule_graph.cpp
===================================================================
--- src/molecule_graph.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule_graph.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -135,8 +135,6 @@
     Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl;
     AtomMap = Calloc<atom *> (AtomCount, "molecule::CreateAdjacencyList - **AtomCount");
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      AtomMap[Walker->nr] = Walker;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+      AtomMap[(*iter)->nr] = (*iter);
     }
 
@@ -496,6 +494,6 @@
   DepthFirstSearchAnalysis_Init(DFS, this);
 
-  DFS.Root = start->next;
-  while (DFS.Root != end) { // if there any atoms at all
+  for (molecule::const_iterator iter = begin(); iter != end();) {
+    DFS.Root = *iter;
     // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all
     DFS.AtomStack->ClearStack();
@@ -537,8 +535,8 @@
 
     // step on to next root
-    while ((DFS.Root != end) && (DFS.Root->GraphNr != -1)) {
-      //Log() << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
-      if (DFS.Root->GraphNr != -1) // if already discovered, step on
-        DFS.Root = DFS.Root->next;
+    while ((iter != end()) && ((*iter)->GraphNr != -1)) {
+      //Log() << Verbose(1) << "Current next subgraph root candidate is " << (*iter)->Name << "." << endl;
+      if ((*iter)->GraphNr != -1) // if already discovered, step on
+        iter++;
     }
   }
@@ -843,7 +841,6 @@
   if (MinRingSize != -1) { // if rings are present
     // go over all atoms
-    Root = mol->start;
-    while (Root->next != mol->end) {
-      Root = Root->next;
+    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      Root = *iter;
 
       if (MinimumRingSize[Root->GetTrueFather()->nr] == mol->AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
@@ -1352,10 +1349,8 @@
   // fill parent list with sons
   Log() << Verbose(3) << "Filling Parent List." << endl;
-  atom *Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    ParentList[Walker->father->nr] = Walker;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    ParentList[(*iter)->father->nr] = (*iter);
     // Outputting List for debugging
-    Log() << Verbose(4) << "Son[" << Walker->father->nr << "] of " << Walker->father << " is " << ParentList[Walker->father->nr] << "." << endl;
+    Log() << Verbose(4) << "Son[" << (*iter)->father->nr << "] of " << (*iter)->father << " is " << ParentList[(*iter)->father->nr] << "." << endl;
   }
 
@@ -1372,20 +1367,17 @@
 {
   bool status = true;
-  atom *Walker = NULL;
   atom *OtherAtom = NULL;
   // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
   Log() << Verbose(3) << "Creating bonds." << endl;
-  Walker = Father->start;
-  while (Walker->next != Father->end) {
-    Walker = Walker->next;
-    if (ParentList[Walker->nr] != NULL) {
-      if (ParentList[Walker->nr]->father != Walker) {
+  for (molecule::const_iterator iter = Father->begin(); iter != Father->end(); ++iter) {
+    if (ParentList[(*iter)->nr] != NULL) {
+      if (ParentList[(*iter)->nr]->father != (*iter)) {
         status = false;
       } else {
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
-          OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
+          OtherAtom = (*Runner)->GetOtherAtom((*iter));
           if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
-            Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
-            mol->AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
+            Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[(*iter)->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
+            mol->AddBond(ParentList[(*iter)->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
           }
         }
Index: src/molecule_pointcloud.cpp
===================================================================
--- src/molecule_pointcloud.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule_pointcloud.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -8,4 +8,5 @@
 #include "atom.hpp"
 #include "config.hpp"
+#include "info.hpp"
 #include "memoryallocator.hpp"
 #include "molecule.hpp"
@@ -31,76 +32,48 @@
 };
 
-/** Return current atom in the list.
- * \return pointer to atom or NULL if none present
+
+/** PointCloud implementation of GoPoint
+ * Uses atoms and STL stuff.
  */
-TesselPoint *molecule::GetPoint() const
+TesselPoint* molecule::GetPoint() const
 {
-  if ((InternalPointer != start) && (InternalPointer != end))
-    return InternalPointer;
-  else
-    return NULL;
+  Info FunctionInfo(__func__);
+  return (*InternalPointer);
 };
 
-/** Return pointer to one after last atom in the list.
- * \return pointer to end marker
- */
-TesselPoint *molecule::GetTerminalPoint() const
-{
-  return end;
-};
-
-/** Return the greatest index of all atoms in the list.
- * \return greatest index
- */
-int molecule::GetMaxId() const
-{
-  return last_atom;
-};
-
-/** Go to next atom.
- * Stops at last one.
+/** PointCloud implementation of GoToNext.
+ * Uses atoms and STL stuff.
  */
 void molecule::GoToNext() const
 {
-  if (InternalPointer != end)
-    InternalPointer = InternalPointer->next;
+  Info FunctionInfo(__func__);
+  if (InternalPointer != end())
+    InternalPointer++;
 };
 
-/** Go to previous atom.
- * Stops at first one.
- */
-void molecule::GoToPrevious() const
-{
-  if (InternalPointer->previous != start)
-    InternalPointer = InternalPointer->previous;
-};
-
-/** Goes to first atom.
+/** PointCloud implementation of GoToFirst.
+ * Uses atoms and STL stuff.
  */
 void molecule::GoToFirst() const
 {
-  InternalPointer = start->next;
+  Info FunctionInfo(__func__);
+  InternalPointer = begin();
 };
 
-/** Goes to last atom.
- */
-void molecule::GoToLast() const
-{
-  InternalPointer = end->previous;
-};
-
-/** Checks whether we have any atoms in molecule.
- * \return true - no atoms, false - not empty
+/** PointCloud implementation of IsEmpty.
+ * Uses atoms and STL stuff.
  */
 bool molecule::IsEmpty() const
 {
-  return (start->next == end);
+  Info FunctionInfo(__func__);
+  return (empty());
 };
 
-/** Checks whether we are at the last atom
- * \return true - current atom is last one, false - is not last one
+/** PointCloud implementation of IsLast.
+ * Uses atoms and STL stuff.
  */
 bool molecule::IsEnd() const
 {
-  return (InternalPointer == end);
+  Info FunctionInfo(__func__);
+  return (InternalPointer == end());
 };
Index: src/molecule_template.hpp
===================================================================
--- src/molecule_template.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/molecule_template.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -24,16 +24,12 @@
 template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() ) const
     {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)();
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)();
   }
 };
 template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() const ) const
     {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)();
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)();
   }
 };
@@ -41,16 +37,12 @@
 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)(t);
   }
 };
 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)(t);
   }
 };
@@ -58,16 +50,12 @@
 template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t, u);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)(t, u);
   }
 };
 template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t, u);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)(t, u);
   }
 };
@@ -75,16 +63,12 @@
 template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t, u, v);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)(t, u, v);
   }
 };
 template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t, u, v);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (((*iter)->node)->*f)(t, u, v);
   }
 };
@@ -96,8 +80,6 @@
 {
   res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)();
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    result += ((*iter)->*f)();
   }
   return result;
@@ -106,8 +88,6 @@
 {
   res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)();
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    result += ((*iter)->*f)();
   }
   return result;
@@ -117,8 +97,6 @@
 {
   res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)(t);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    result += ((*iter)->*f)(t);
   }
   return result;
@@ -127,8 +105,6 @@
 {
   res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)(t);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    result += ((*iter)->*f)(t);
   }
   return result;
@@ -141,16 +117,12 @@
 template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *)) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*f)((*iter));
   }
 };
 template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *) const) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*f)((*iter));
   }
 };
@@ -161,16 +133,12 @@
 template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (copy->*f)((*iter));
   }
 };
 template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (copy->*f)((*iter));
   }
 };
@@ -181,36 +149,28 @@
 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)())
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)())
+      (copy->*f)((*iter));
   }
 };
 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () const ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)())
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)())
+      (copy->*f)((*iter));
   }
 };
 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)())
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)())
+      (copy->*f)((*iter));
   }
 };
 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) () const ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)())
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)())
+      (copy->*f)((*iter));
   }
 };
@@ -218,36 +178,28 @@
 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T) const, T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T), T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T) const, T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t))
+      (copy->*f)((*iter));
   }
 };
@@ -255,36 +207,28 @@
 template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u))
+      (copy->*f)((*iter));
   }
 };
@@ -292,36 +236,28 @@
 template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u,v))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u,v))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u,v))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u,v))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u,v))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u,v))
+      (copy->*f)((*iter));
   }
 };
 template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u,v))
-      (copy->*f)(Walker);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    if (((*iter)->*condition)(t,u,v))
+      (copy->*f)((*iter));
   }
 };
@@ -332,16 +268,12 @@
 template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)()) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)();
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)();
   }
 };
 template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)();
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)();
   }
 };
@@ -349,16 +281,12 @@
 template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t);
   }
 };
 template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t);
   }
 };
@@ -366,16 +294,12 @@
 template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U), T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t, u);
   }
 };
 template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t, u);
   }
 };
@@ -383,16 +307,12 @@
 template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t, u, v);
   }
 };
 template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t, u, v);
   }
 };
@@ -400,16 +320,12 @@
 template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v, w);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t, u, v, w);
   }
 };
 template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v, w);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    ((*iter)->*f)(t, u, v, w);
   }
 };
@@ -420,25 +336,19 @@
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *) ) const
 {
-  atom *Walker = start;
   int inc = 1;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Setor) (&array[(Walker->*index)], &inc);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*Setor) (&array[((*iter)->*index)], &inc);
   }
 };
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Setor) (&array[(Walker->*index)], &value);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*Setor) (&array[((*iter)->*index)], &value);
   }
 };
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T *value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Setor) (&array[(Walker->*index)], value);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*Setor) (&array[((*iter)->*index)], value);
   }
 };
@@ -446,25 +356,19 @@
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *) ) const
 {
-  atom *Walker = start;
   int inc = 1;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Setor) (&array[(Walker->type->*index)], &inc);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*Setor) (&array[((*iter)->type->*index)], &inc);
   }
 };
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Setor) (&array[(Walker->type->*index)], &value);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*Setor) (&array[((*iter)->type->*index)], &value);
   }
 };
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T *value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Setor) (&array[(Walker->type->*index)], value);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*Setor) (&array[((*iter)->type->*index)], value);
   }
 };
@@ -472,41 +376,31 @@
 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ atom::*value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    array[(Walker->*index)] = (Walker->*Setor) (Walker->*value);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    array[((*iter)->*index)] = ((*iter)->*Setor) ((*iter)->*value);
   }
 };
 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ atom::*value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    array[(Walker->*index)] = (Walker->*Setor) (Walker->*value);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    array[((*iter)->*index)] = ((*iter)->*Setor) ((*iter)->*value);
   }
 };
 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ &vect ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    array[(Walker->*index)] = (Walker->*Setor) (vect);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    array[((*iter)->*index)] = ((*iter)->*Setor) (vect);
   }
 };
 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ &vect ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    array[(Walker->*index)] = (Walker->*Setor) (vect);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    array[((*iter)->*index)] = ((*iter)->*Setor) (vect);
   }
 };
 template <typename T, typename typ, typename typ2> void molecule::SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    Walker->*value = array[(Walker->*index)];
-    //Log() << Verbose(2) << *Walker << " gets " << (Walker->*value); << endl;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*iter)->*value = array[((*iter)->*index)];
+    //Log() << Verbose(2) << *(*iter) << " gets " << ((*iter)->*value); << endl;
   }
 };
@@ -514,9 +408,7 @@
 template <typename T, typename typ> void molecule::SetAtomValueToValue ( T value, T typ::*ptr ) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    Walker->*ptr = value;
-    //Log() << Verbose(2) << *Walker << " gets " << (Walker->*ptr) << endl;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*iter)->*ptr = value;
+    //Log() << Verbose(2) << *(*iter) << " gets " << ((*iter)->*ptr) << endl;
   }
 };
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/moleculelist.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -68,6 +68,4 @@
   int Count, Counter, aCounter, bCounter;
   int flag;
-  atom *aWalker = NULL;
-  atom *bWalker = NULL;
 
   // sort each atom list and put the numbers into a list, then go through
@@ -84,28 +82,27 @@
 
       // fill the lists
-      aWalker = (**(molecule **) a).start;
-      bWalker = (**(molecule **) b).start;
       Counter = 0;
       aCounter = 0;
       bCounter = 0;
-      while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
-        aWalker = aWalker->next;
-        bWalker = bWalker->next;
-        if (aWalker->GetTrueFather() == NULL)
+      molecule::const_iterator aiter = (**(molecule **) a).begin();
+      molecule::const_iterator biter = (**(molecule **) b).begin();
+      for (;(aiter != (**(molecule **) a).end()) && (biter != (**(molecule **) b).end());
+          ++aiter, ++biter) {
+        if ((*aiter)->GetTrueFather() == NULL)
           aList[Counter] = Count + (aCounter++);
         else
-          aList[Counter] = aWalker->GetTrueFather()->nr;
-        if (bWalker->GetTrueFather() == NULL)
+          aList[Counter] = (*aiter)->GetTrueFather()->nr;
+        if ((*biter)->GetTrueFather() == NULL)
           bList[Counter] = Count + (bCounter++);
         else
-          bList[Counter] = bWalker->GetTrueFather()->nr;
+          bList[Counter] = (*biter)->GetTrueFather()->nr;
         Counter++;
       }
       // check if AtomCount was for real
       flag = 0;
-      if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
+      if ((aiter == (**(molecule **) a).end()) && (biter != (**(molecule **) b).end())) {
         flag = -1;
       } else {
-        if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end))
+        if ((aiter != (**(molecule **) a).end()) && (biter == (**(molecule **) b).end()))
           flag = 1;
       }
@@ -142,5 +139,4 @@
 {
   element* Elemental = NULL;
-  atom *Walker = NULL;
   int Counts[MAX_ELEMENTS];
   double size=0;
@@ -160,10 +156,8 @@
       // count atoms per element and determine size of bounding sphere
       size=0.;
-      Walker = (*ListRunner)->start;
-      while (Walker->next != (*ListRunner)->end) {
-        Walker = Walker->next;
-        Counts[Walker->type->Z]++;
-        if (Walker->x.DistanceSquared(&Origin) > size)
-          size = Walker->x.DistanceSquared(&Origin);
+      for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
+        Counts[(*iter)->type->Z]++;
+        if ((*iter)->x.DistanceSquared(&Origin) > size)
+          size = (*iter)->x.DistanceSquared(&Origin);
       }
       // output Index, Name, number of atoms, chemical formula
@@ -204,11 +198,9 @@
 
   // put all molecules of src into mol
-  atom *Walker = srcmol->start;
-  atom *NextAtom = Walker->next;
-  while (NextAtom != srcmol->end) {
-    Walker = NextAtom;
-    NextAtom = Walker->next;
-    srcmol->UnlinkAtom(Walker);
-    mol->AddAtom(Walker);
+  molecule::iterator runner;
+  for (molecule::iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
+    runner = iter++;
+    srcmol->UnlinkAtom((*runner));
+    mol->AddAtom((*runner));
   }
 
@@ -230,10 +222,7 @@
 
   // put all molecules of src into mol
-  atom *Walker = srcmol->start;
-  atom *NextAtom = Walker->next;
-  while (NextAtom != srcmol->end) {
-    Walker = NextAtom;
-    NextAtom = Walker->next;
-    Walker = mol->AddCopyAtom(Walker);
+  atom *Walker = NULL;
+  for (molecule::iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
+    Walker = mol->AddCopyAtom((*iter));
     Walker->father = Walker;
   }
@@ -338,12 +327,10 @@
 
   // for each of the source atoms check whether we are in- or outside and add copy atom
-  atom *Walker = srcmol->start;
   int nr=0;
-  while (Walker->next != srcmol->end) {
-    Walker = Walker->next;
-    Log() << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl;
-    if (!TesselStruct->IsInnerPoint(Walker->x, LCList)) {
-      CopyAtoms[Walker->nr] = Walker->clone();
-      mol->AddAtom(CopyAtoms[Walker->nr]);
+  for (molecule::const_iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
+    Log() << Verbose(2) << "INFO: Current Walker is " << *(*iter) << "." << endl;
+    if (!TesselStruct->IsInnerPoint((*iter)->x, LCList)) {
+      CopyAtoms[(*iter)->nr] = (*iter)->clone();
+      mol->AddAtom(CopyAtoms[(*iter)->nr]);
       nr++;
     } else {
@@ -384,6 +371,4 @@
 bool MoleculeListClass::AddHydrogenCorrection(char *path)
 {
-  atom *Walker = NULL;
-  atom *Runner = NULL;
   bond *Binder = NULL;
   double ***FitConstant = NULL, **correction = NULL;
@@ -490,20 +475,16 @@
         correction[k][j] = 0.;
     // 2. take every hydrogen that is a saturated one
-    Walker = (*ListRunner)->start;
-    while (Walker->next != (*ListRunner)->end) {
-      Walker = Walker->next;
-      //Log() << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
-      if ((Walker->type->Z == 1) && ((Walker->father == NULL)
-          || (Walker->father->type->Z != 1))) { // if it's a hydrogen
-        Runner = (*ListRunner)->start;
-        while (Runner->next != (*ListRunner)->end) {
-          Runner = Runner->next;
-          //Log() << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
+    for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
+      //Log() << Verbose(1) << "(*iter): " << *(*iter) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
+      if (((*iter)->type->Z == 1) && (((*iter)->father == NULL)
+          || ((*iter)->father->type->Z != 1))) { // if it's a hydrogen
+        for (molecule::const_iterator runner = (*ListRunner)->begin(); runner != (*ListRunner)->end(); ++runner) {
+          //Log() << Verbose(2) << "Runner: " << *(*runner) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
-          Binder = *(Runner->ListOfBonds.begin());
-          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (Binder->GetOtherAtom(Runner) != Binder->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
+          Binder = *((*runner)->ListOfBonds.begin());
+          if (((*runner)->type->Z == 1) && ((*runner)->nr > (*iter)->nr) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
             // 4. evaluate the morse potential for each matrix component and add up
-            distance = Runner->x.Distance(&Walker->x);
-            //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
+            distance = (*runner)->x.Distance(&(*iter)->x);
+            //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *(*runner) << "<= " << distance << "=>" << *(*iter) << ":" << endl;
             for (int k = 0; k < a; k++) {
               for (int j = 0; j < b; j++) {
@@ -579,5 +560,4 @@
   ofstream ForcesFile;
   stringstream line;
-  atom *Walker = NULL;
   element *runner = NULL;
 
@@ -594,11 +574,9 @@
         runner = runner->next;
         if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms
-          Walker = (*ListRunner)->start;
-          while (Walker->next != (*ListRunner)->end) { // go through every atom of this element
-            Walker = Walker->next;
-            if (Walker->type->Z == runner->Z) {
-              if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
-                //Log() << Verbose(0) << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
-                ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t";
+          for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
+            if ((*iter)->type->Z == runner->Z) {
+              if (((*iter)->GetTrueFather() != NULL) && ((*iter)->GetTrueFather() != (*iter))) {// if there is a rea
+                //Log() << Verbose(0) << "Walker is " << *(*iter) << " with true father " << *( (*iter)->GetTrueFather()) << ", it
+                ForcesFile << SortIndex[(*iter)->GetTrueFather()->nr] << "\t";
               } else
                 // otherwise a -1 to indicate an added saturation hydrogen
@@ -636,5 +614,4 @@
   bool result = true;
   bool intermediateResult = true;
-  atom *Walker = NULL;
   Vector BoxDimension;
   char *FragmentNumber = NULL;
@@ -672,8 +649,6 @@
     // list atoms in fragment for debugging
     Log() << Verbose(2) << "Contained atoms: ";
-    Walker = (*ListRunner)->start;
-    while (Walker->next != (*ListRunner)->end) {
-      Walker = Walker->next;
-      Log() << Verbose(0) << Walker->Name << " ";
+    for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
+      Log() << Verbose(0) << (*iter)->Name << " ";
     }
     Log() << Verbose(0) << endl;
@@ -751,6 +726,4 @@
 {
   molecule *mol = World::get()->createMolecule();
-  atom *Walker = NULL;
-  atom *Advancer = NULL;
   bond *Binder = NULL;
   bond *Stepper = NULL;
@@ -758,12 +731,9 @@
   for (MoleculeList::iterator MolRunner = ListOfMolecules.begin(); !ListOfMolecules.empty(); MolRunner = ListOfMolecules.begin()) {
     // shift all atoms to new molecule
-    Advancer = (*MolRunner)->start->next;
-    while (Advancer != (*MolRunner)->end) {
-      Walker = Advancer;
-      Advancer = Advancer->next;
-      Log() << Verbose(3) << "Re-linking " << *Walker << "..." << endl;
-      unlink(Walker);
-      Walker->father = Walker;
-      mol->AddAtom(Walker);    // counting starts at 1
+    for (molecule::iterator iter = (*MolRunner)->begin(); (*MolRunner)->empty(); iter = (*MolRunner)->begin()) {
+      Log() << Verbose(3) << "Re-linking " << *(*iter) << "..." << endl;
+      (*iter)->father = (*iter);
+      mol->AddAtom((*iter));    // counting starts at 1
+      (*MolRunner)->erase(iter);
     }
     // remove all bonds
@@ -812,11 +782,8 @@
   int *MolMap = Calloc<int>(mol->AtomCount, "config::Load() - *MolMap");
   MoleculeLeafClass *MolecularWalker = Subgraphs;
-  Walker = NULL;
   while (MolecularWalker->next != NULL) {
     MolecularWalker = MolecularWalker->next;
-    Walker = MolecularWalker->Leaf->start;
-    while (Walker->next != MolecularWalker->Leaf->end) {
-      Walker = Walker->next;
-      MolMap[Walker->GetTrueFather()->nr] = FragmentCounter+1;
+    for (molecule::const_iterator iter = MolecularWalker->Leaf->begin(); iter != MolecularWalker->Leaf->end(); ++iter) {
+      MolMap[(*iter)->GetTrueFather()->nr] = FragmentCounter+1;
     }
     FragmentCounter++;
@@ -824,18 +791,16 @@
 
   // 4c. relocate atoms to new molecules and remove from Leafs
-  Walker = NULL;
-  while (mol->start->next != mol->end) {
-    Walker = mol->start->next;
-    if ((Walker->nr <0) || (Walker->nr >= mol->AtomCount)) {
-      eLog() << Verbose(0) << "Index of atom " << *Walker << " is invalid!" << endl;
+  for (molecule::iterator iter = mol->begin(); !mol->empty(); iter = mol->begin()) {
+    if (((*iter)->nr <0) || ((*iter)->nr >= mol->AtomCount)) {
+      eLog() << Verbose(0) << "Index of atom " << *(*iter) << " is invalid!" << endl;
       performCriticalExit();
     }
-    FragmentCounter = MolMap[Walker->nr];
+    FragmentCounter = MolMap[(*iter)->nr];
     if (FragmentCounter != 0) {
-      Log() << Verbose(3) << "Re-linking " << *Walker << "..." << endl;
-      unlink(Walker);
-      molecules[FragmentCounter-1]->AddAtom(Walker);    // counting starts at 1
+      Log() << Verbose(3) << "Re-linking " << *(*iter) << "..." << endl;
+      molecules[FragmentCounter-1]->AddAtom((*iter));    // counting starts at 1
+      mol->erase(iter);
     } else {
-      eLog() << Verbose(0) << "Atom " << *Walker << " not associated to molecule!" << endl;
+      eLog() << Verbose(0) << "Atom " << *(*iter) << " not associated to molecule!" << endl;
       performCriticalExit();
     }
@@ -845,5 +810,5 @@
   while (mol->first->next != mol->last) {
     Binder = mol->first->next;
-    Walker = Binder->leftatom;
+    const atom * const Walker = Binder->leftatom;
     unlink(Binder);
     link(Binder,molecules[MolMap[Walker->nr]-1]->last);   // counting starts at 1
@@ -867,12 +832,7 @@
 int MoleculeListClass::CountAllAtoms() const
 {
-  atom *Walker = NULL;
   int AtomNo = 0;
   for (MoleculeList::const_iterator MolWalker = ListOfMolecules.begin(); MolWalker != ListOfMolecules.end(); MolWalker++) {
-    Walker = (*MolWalker)->start;
-    while (Walker->next != (*MolWalker)->end) {
-      Walker = Walker->next;
-      AtomNo++;
-    }
+    AtomNo += (*MolWalker)->size();
   }
   return AtomNo;
@@ -1049,5 +1009,4 @@
 bool MoleculeLeafClass::FillBondStructureFromReference(const molecule * const reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
 {
-  atom *Walker = NULL;
   atom *OtherWalker = NULL;
   atom *Father = NULL;
@@ -1073,16 +1032,14 @@
     }
 
-    Walker = Leaf->start;
-    while (Walker->next != Leaf->end) {
-      Walker = Walker->next;
-      Father = Walker->GetTrueFather();
+    for(molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ++iter) {
+      Father = (*iter)->GetTrueFather();
       AtomNo = Father->nr; // global id of the current walker
       for (BondList::const_iterator Runner = Father->ListOfBonds.begin(); Runner != Father->ListOfBonds.end(); (++Runner)) {
-        OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
+        OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr]; // local copy of current bond partner of walker
         if (OtherWalker != NULL) {
-          if (OtherWalker->nr > Walker->nr)
-            Leaf->AddBond(Walker, OtherWalker, (*Runner)->BondDegree);
+          if (OtherWalker->nr > (*iter)->nr)
+            Leaf->AddBond((*iter), OtherWalker, (*Runner)->BondDegree);
         } else {
-          Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
+          Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr << "] is NULL!" << endl;
           status = false;
         }
@@ -1111,5 +1068,5 @@
 bool MoleculeLeafClass::FillRootStackForSubgraphs(KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
 {
-  atom *Walker = NULL, *Father = NULL;
+  atom *Father = NULL;
 
   if (RootStack != NULL) {
@@ -1117,13 +1074,11 @@
     if (&(RootStack[FragmentCounter]) != NULL) {
       RootStack[FragmentCounter].clear();
-      Walker = Leaf->start;
-      while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
-        Walker = Walker->next;
-        Father = Walker->GetTrueFather();
+      for(molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ++iter) {
+        Father = (*iter)->GetTrueFather();
         if (AtomMask[Father->nr]) // apply mask
 #ifdef ADDHYDROGEN
-          if (Walker->type->Z != 1) // skip hydrogen
+          if ((*iter)->type->Z != 1) // skip hydrogen
 #endif
-          RootStack[FragmentCounter].push_front(Walker->nr);
+          RootStack[FragmentCounter].push_front((*iter)->nr);
       }
       if (next != NULL)
@@ -1164,5 +1119,5 @@
 
   if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
-    status = status && CreateFatherLookupTable(Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
+    status = status && Leaf->CreateFatherLookupTable(ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
     FreeList = FreeList && true;
   }
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/tesselation.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -1100,15 +1100,4 @@
 };
 
-/** PointCloud implementation of GetTerminalPoint.
- * Uses PointsOnBoundary and STL stuff.
- */   
-TesselPoint * Tesselation::GetTerminalPoint() const
-{
-	Info FunctionInfo(__func__);
-  PointMap::const_iterator Runner = PointsOnBoundary.end();
-  Runner--;
-  return (Runner->second->node);
-};
-
 /** PointCloud implementation of GoToNext.
  * Uses PointsOnBoundary and STL stuff.
@@ -1121,14 +1110,4 @@
 };
 
-/** PointCloud implementation of GoToPrevious.
- * Uses PointsOnBoundary and STL stuff.
- */   
-void Tesselation::GoToPrevious() const
-{
-	Info FunctionInfo(__func__);
-  if (InternalPointer != PointsOnBoundary.begin())
-    InternalPointer--;
-};
-
 /** PointCloud implementation of GoToFirst.
  * Uses PointsOnBoundary and STL stuff.
@@ -1138,14 +1117,4 @@
 	Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.begin();
-};
-
-/** PointCloud implementation of GoToLast.
- * Uses PointsOnBoundary and STL stuff.
- */
-void Tesselation::GoToLast() const
-{
-	Info FunctionInfo(__func__);
-  InternalPointer = PointsOnBoundary.end();
-  InternalPointer--;
 };
 
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/tesselation.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -234,10 +234,7 @@
   virtual Vector *GetCenter() const { return NULL; };
   virtual TesselPoint *GetPoint() const { return NULL; };
-  virtual TesselPoint *GetTerminalPoint() const { return NULL; };
   virtual int GetMaxId() const { return 0; };
   virtual void GoToNext() const {};
-  virtual void GoToPrevious() const {};
   virtual void GoToFirst() const {};
-  virtual void GoToLast() const {};
   virtual bool IsEmpty() const { return true; };
   virtual bool IsEnd() const { return true; };
@@ -342,9 +339,6 @@
     virtual Vector *GetCenter(ofstream *out) const;
     virtual TesselPoint *GetPoint() const;
-    virtual TesselPoint *GetTerminalPoint() const;
     virtual void GoToNext() const;
-    virtual void GoToPrevious() const;
     virtual void GoToFirst() const;
-    virtual void GoToLast() const;
     virtual bool IsEmpty() const;
     virtual bool IsEnd() const;
Index: src/unittests/AnalysisCorrelationToPointUnitTest.cpp
===================================================================
--- src/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -25,5 +25,4 @@
 #include "periodentafel.hpp"
 #include "tesselation.hpp"
-#include "World.hpp"
 
 #ifdef HAVE_TESTRUNNER
Index: src/unittests/analysisbondsunittest.cpp
===================================================================
--- src/unittests/analysisbondsunittest.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/unittests/analysisbondsunittest.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -25,4 +25,5 @@
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "World.hpp"
 
 #ifdef HAVE_TESTRUNNER
Index: src/unittests/bondgraphunittest.cpp
===================================================================
--- src/unittests/bondgraphunittest.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/unittests/bondgraphunittest.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -110,4 +110,12 @@
 };
 
+/** Tests whether setup worked.
+ */
+void BondGraphTest::SetupTest()
+{
+  CPPUNIT_ASSERT_EQUAL (false, TestMolecule->empty());
+  CPPUNIT_ASSERT_EQUAL ((size_t)4, TestMolecule->size());
+};
+
 /** UnitTest for BondGraphTest::LoadBondLengthTable().
  */
@@ -124,9 +132,9 @@
 void BondGraphTest::ConstructGraphTest()
 {
-  atom *Walker = TestMolecule->start->next;
-  atom *Runner = TestMolecule->end->previous;
-  CPPUNIT_ASSERT( TestMolecule->end != Walker );
+  molecule::iterator Walker = TestMolecule->begin();
+  molecule::iterator Runner = TestMolecule->begin();
+  Runner++;
   CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(*filename) );
   CPPUNIT_ASSERT_EQUAL( true , BG->ConstructBondGraph(TestMolecule) );
-  CPPUNIT_ASSERT_EQUAL( true , Walker->IsBondedTo(Runner) );
+  CPPUNIT_ASSERT_EQUAL( true , (*Walker)->IsBondedTo((*Runner)) );
 };
Index: src/unittests/bondgraphunittest.hpp
===================================================================
--- src/unittests/bondgraphunittest.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/unittests/bondgraphunittest.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -22,4 +22,5 @@
 {
     CPPUNIT_TEST_SUITE( BondGraphTest) ;
+    CPPUNIT_TEST ( SetupTest );
     CPPUNIT_TEST ( LoadTableTest );
     CPPUNIT_TEST ( ConstructGraphTest );
@@ -29,4 +30,5 @@
       void setUp();
       void tearDown();
+      void SetupTest();
       void LoadTableTest();
       void ConstructGraphTest();
Index: src/unittests/listofbondsunittest.cpp
===================================================================
--- src/unittests/listofbondsunittest.cpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/unittests/listofbondsunittest.cpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -90,4 +90,13 @@
 };
 
+/** Tests whether setup worked correctly.
+ *
+ */
+void ListOfBondsTest::SetupTest()
+{
+  CPPUNIT_ASSERT_EQUAL( false, TestMolecule->empty() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, TestMolecule->size() );
+};
+
 /** Unit Test of molecule::AddBond()
  *
@@ -96,6 +105,8 @@
 {
   bond *Binder = NULL;
-  atom *atom1 = TestMolecule->start->next;
-  atom *atom2 = atom1->next;
+  molecule::iterator iter = TestMolecule->begin();
+  atom *atom1 = *iter;
+  iter++;
+  atom *atom2 = *iter;
   CPPUNIT_ASSERT( atom1 != NULL );
   CPPUNIT_ASSERT( atom2 != NULL );
@@ -124,6 +135,8 @@
 {
   bond *Binder = NULL;
-  atom *atom1 = TestMolecule->start->next;
-  atom *atom2 = atom1->next;
+  molecule::iterator iter = TestMolecule->begin();
+  atom *atom1 = *iter;
+  iter++;
+  atom *atom2 = *iter;
   CPPUNIT_ASSERT( atom1 != NULL );
   CPPUNIT_ASSERT( atom2 != NULL );
@@ -150,7 +163,10 @@
 {
   bond *Binder = NULL;
-  atom *atom1 = TestMolecule->start->next;
-  atom *atom2 = atom1->next;
-  atom *atom3 = atom2->next;
+  molecule::iterator iter = TestMolecule->begin();
+  atom *atom1 = *iter;
+  iter++;
+  atom *atom2 = *iter;
+  iter++;
+  atom *atom3 = *iter;
   CPPUNIT_ASSERT( atom1 != NULL );
   CPPUNIT_ASSERT( atom2 != NULL );
@@ -189,6 +205,8 @@
 {
   bond *Binder = NULL;
-  atom *atom1 = TestMolecule->start->next;
-  atom *atom2 = atom1->next;
+  molecule::iterator iter = TestMolecule->begin();
+  atom *atom1 = *iter;
+  iter++;
+  atom *atom2 = *iter;
   CPPUNIT_ASSERT( atom1 != NULL );
   CPPUNIT_ASSERT( atom2 != NULL );
@@ -215,6 +233,8 @@
 {
   bond *Binder = NULL;
-  atom *atom1 = TestMolecule->start->next;
-  atom *atom2 = atom1->next;
+  molecule::iterator iter = TestMolecule->begin();
+  atom *atom1 = *iter;
+  iter++;
+  atom *atom2 = *iter;
   CPPUNIT_ASSERT( atom1 != NULL );
   CPPUNIT_ASSERT( atom2 != NULL );
@@ -240,6 +260,8 @@
 {
   bond *Binder = NULL;
-  atom *atom1 = TestMolecule->start->next;
-  atom *atom2 = atom1->next;
+  molecule::iterator iter = TestMolecule->begin();
+  atom *atom1 = *iter;
+  iter++;
+  atom *atom2 = *iter;
   CPPUNIT_ASSERT( atom1 != NULL );
   CPPUNIT_ASSERT( atom2 != NULL );
Index: src/unittests/listofbondsunittest.hpp
===================================================================
--- src/unittests/listofbondsunittest.hpp	(revision e87acf50f98112d75e2652d26e782186d90e5f4e)
+++ src/unittests/listofbondsunittest.hpp	(revision 9879f6f8e5ad9364ccd32e2d6735f2ffc89cf18d)
@@ -20,4 +20,5 @@
 {
     CPPUNIT_TEST_SUITE( ListOfBondsTest) ;
+    CPPUNIT_TEST ( SetupTest );
     CPPUNIT_TEST ( AddingBondTest );
     CPPUNIT_TEST ( RemovingBondTest );
@@ -31,4 +32,5 @@
       void setUp();
       void tearDown();
+      void SetupTest();
       void AddingBondTest();
       void RemovingBondTest();
