Index: src/Actions/FragmentationAction/ConstructBondGraphAction.cpp
===================================================================
--- src/Actions/FragmentationAction/ConstructBondGraphAction.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/Actions/FragmentationAction/ConstructBondGraphAction.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -32,4 +32,5 @@
 #include "molecule.hpp"
 #include "World.hpp"
+#include "WorldTime.hpp"
 
 #include <iostream>
@@ -67,7 +68,12 @@
   // remove every bond from the selected atoms' list
   int AtomCount = 0;
-  for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection(); AtomRunner != World::getInstance().endAtomSelection(); ++AtomRunner) {
+  for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
+      AtomRunner != World::getInstance().endAtomSelection();
+      ++AtomRunner) {
     AtomCount++;
-    for(BondList::iterator BondRunner = (AtomRunner->second)->ListOfBonds.begin(); !(AtomRunner->second)->ListOfBonds.empty(); BondRunner = (AtomRunner->second)->ListOfBonds.begin())
+    BondList& ListOfBonds = (AtomRunner->second)->getListOfBonds();
+    for(BondList::iterator BondRunner = ListOfBonds.begin();
+        !ListOfBonds.empty();
+        BondRunner = ListOfBonds.begin())
       if ((*BondRunner)->leftatom == AtomRunner->second)
         delete((*BondRunner));
@@ -81,5 +87,7 @@
     DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
     LinkedCell::LinkedNodes list;
-    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection(); AtomRunner != World::getInstance().endAtomSelection(); ++AtomRunner) {
+    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
+        AtomRunner != World::getInstance().endAtomSelection();
+        ++AtomRunner) {
       list.push_back(AtomRunner->second);
     }
@@ -92,5 +100,7 @@
     std::map<TesselPoint *, int> AtomIds;
     int i=0;
-    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection(); AtomRunner != World::getInstance().endAtomSelection(); ++AtomRunner) {
+    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
+        AtomRunner != World::getInstance().endAtomSelection();
+        ++AtomRunner) {
       AtomIds.insert(pair<TesselPoint *, int> (AtomRunner->second, i++) );
     }
@@ -104,5 +114,7 @@
 //          Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
           if (List != NULL) {
-            for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin();
+                Runner != List->end();
+                Runner++) {
               Walker = dynamic_cast<atom*>(*Runner);
               ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
Index: src/Actions/FragmentationAction/SubgraphDissectionAction.cpp
===================================================================
--- src/Actions/FragmentationAction/SubgraphDissectionAction.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/Actions/FragmentationAction/SubgraphDissectionAction.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -78,7 +78,11 @@
     vector <atom *> allatoms = World::getInstance().getAllAtoms();
     for(vector<atom *>::iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) {
-      BondCount += (*AtomRunner)->ListOfBonds.size();
-//      for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
+      const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+      BondCount += ListOfBonds.size();
+//      for(BondList::iterator BondRunner = ListOfBonds.begin();
+//          !ListOfBonds.empty();
+//          BondRunner = ListOfBonds.begin()) {
 //        delete(*BondRunner);
+//      }
       mol->AddAtom(*AtomRunner);
     }
Index: src/Actions/MoleculeAction/CopyAction.cpp
===================================================================
--- src/Actions/MoleculeAction/CopyAction.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/Actions/MoleculeAction/CopyAction.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -59,7 +59,13 @@
   MoleculeCopyState *state = assert_cast<MoleculeCopyState*>(_state.get());
 
-  for (molecule::iterator AtomRunner = state->copy->begin(); !state->copy->empty(); AtomRunner = state->copy->begin()) {
-    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
+  for (molecule::iterator AtomRunner = state->copy->begin();
+      !state->copy->empty();
+      AtomRunner = state->copy->begin()) {
+    BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::iterator BondRunner = ListOfBonds.begin();
+        !ListOfBonds.empty();
+        BondRunner = ListOfBonds.begin()) {
       delete(*BondRunner);
+    }
     atom *Walker = *AtomRunner;
     state->copy->erase(AtomRunner);
Index: src/Parser/PdbParser.cpp
===================================================================
--- src/Parser/PdbParser.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/Parser/PdbParser.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -325,4 +325,6 @@
  * Writes the neighbor information of one atom to the provided stream.
  *
+ * Note that ListOfBonds of WorldTime::CurrentTime is used.
+ *
  * \param *file  where to write neighbor information to
  * \param MaxnumberOfNeighbors of neighbors
@@ -331,6 +333,7 @@
 void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) {
   int MaxNo = MaxnumberOfNeighbors;
-  if (!currentAtom->ListOfBonds.empty()) {
-    for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) {
+  const BondList & ListOfBonds = currentAtom->getListOfBonds();
+  if (!ListOfBonds.empty()) {
+    for(BondList::const_iterator currentBond = ListOfBonds.begin(); currentBond != ListOfBonds.end(); ++currentBond) {
       if (MaxNo >= MaxnumberOfNeighbors) {
         *file << "CONECT";
Index: src/Parser/TremoloParser.cpp
===================================================================
--- src/Parser/TremoloParser.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/Parser/TremoloParser.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -26,4 +26,5 @@
 #include "TremoloParser.hpp"
 #include "World.hpp"
+#include "WorldTime.hpp"
 #include "atom.hpp"
 #include "bond.hpp"
@@ -258,4 +259,6 @@
  * Writes the neighbor information of one atom to the provided stream.
  *
+ * Note that ListOfBonds of WorldTime::CurrentTime is used.
+ *
  * \param stream where to write neighbor information to
  * \param number of neighbors
@@ -263,9 +266,10 @@
  */
 void TremoloParser::writeNeighbors(ostream* file, int numberOfNeighbors, atom* currentAtom) {
-  BondList::iterator currentBond = currentAtom->ListOfBonds.begin();
+  const BondList& ListOfBonds = currentAtom->getListOfBonds();
+  BondList::const_iterator currentBond = ListOfBonds.begin();
   for (int i = 0; i < numberOfNeighbors; i++) {
-    *file << (currentBond != currentAtom->ListOfBonds.end()
+    *file << (currentBond != ListOfBonds.end()
         ? (*currentBond)->GetOtherAtom(currentAtom)->getId()+1 : 0) << "\t";
-    if (currentBond != currentAtom->ListOfBonds.end())
+    if (currentBond != ListOfBonds.end())
       ++currentBond;
   }
Index: src/WorldTime.hpp
===================================================================
--- src/WorldTime.hpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/WorldTime.hpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -22,4 +22,6 @@
   //!> own unit test should be friend to access private variables.
   friend class WorldTimeTest;
+  //!> ListOfBonds test has to be friend to access private variables.
+  friend class ListOfBondsTest;
 
   //!> World is friend to access setTime (alone!)
Index: src/analysis_bonds.cpp
===================================================================
--- src/analysis_bonds.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/analysis_bonds.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -43,5 +43,6 @@
   int AtomCount = 0;
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
-    const int count = (*iter)->ListOfBonds.size();
+    const BondList& ListOfBonds = (*iter)->getListOfBonds();
+    const int count = ListOfBonds.size();
     if (Max < count)
       Max = count;
@@ -72,6 +73,9 @@
   int AtomNo = 0;
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
-    if ((*iter)->getType() == type1)
-      for (BondList::const_iterator BondRunner = (*iter)->ListOfBonds.begin(); BondRunner != (*iter)->ListOfBonds.end(); BondRunner++)
+    if ((*iter)->getType() == type1) {
+      const BondList& ListOfBonds = (*iter)->getListOfBonds();
+      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+          BondRunner != ListOfBonds.end();
+          BondRunner++)
         if ((*BondRunner)->GetOtherAtom((*iter))->getType() == type2) {
           const double distance = (*BondRunner)->GetDistanceSquared();
@@ -83,4 +87,5 @@
           AtomNo++;
         }
+    }
   }
   if (Max < 0) {
@@ -161,5 +166,8 @@
               InterfaceFlag = (InterfaceElement == NULL);
               Interface2Flag = (Interface2Element == NULL);
-              for (BondList::const_iterator BondRunner = (*Runner)->ListOfBonds.begin(); BondRunner != (*Runner)->ListOfBonds.end(); BondRunner++) {
+              const BondList& ListOfBonds = (*Runner)->getListOfBonds();
+              for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+                  BondRunner != ListOfBonds.end();
+                  BondRunner++) {
                 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
                 // if hydrogen, check angle to be greater(!) than 30 degrees
@@ -187,5 +195,8 @@
               if (InterfaceFlag && Interface2Flag && OtherHydrogenFlag) {
                 // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
-                for (BondList::const_iterator BondRunner = (*Walker)->ListOfBonds.begin(); BondRunner != (*Walker)->ListOfBonds.end(); BondRunner++) {
+                const BondList& ListOfBonds = (*Walker)->getListOfBonds();
+                for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+                    BondRunner != ListOfBonds.end();
+                    BondRunner++) {
                   atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
                   if (OtherAtom->getType()->getAtomicNumber() == 1) {
@@ -223,5 +234,8 @@
       atom * theAtom = *Walker;
       if ((theAtom->getType() == first) || (theAtom->getType() == second)) {  // first element matches
-        for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
+        const BondList& ListOfBonds = theAtom->getListOfBonds();
+        for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+            BondRunner != ListOfBonds.end();
+            BondRunner++) {
           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
           if (((OtherAtom->getType() == first) || (OtherAtom->getType() == second)) && (theAtom->nr < OtherAtom->nr)) {
@@ -261,5 +275,8 @@
         for (int i=0;i<2;i++)
           MatchFlag[i] = false;
-        for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
+        const BondList& ListOfBonds = theAtom->getListOfBonds();
+        for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+            BondRunner != ListOfBonds.end();
+            BondRunner++) {
           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
           for (int i=0;i<2;i++)
Index: src/analysis_correlation.cpp
===================================================================
--- src/analysis_correlation.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/analysis_correlation.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -66,6 +66,7 @@
       ++atomiter) {
     // go through all bonds
-    for (BondList::const_iterator bonditer = (*atomiter)->ListOfBonds.begin();
-        bonditer != (*atomiter)->ListOfBonds.end();
+    const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
+    for (BondList::const_iterator bonditer = ListOfBonds.begin();
+        bonditer != ListOfBonds.end();
         ++bonditer) {
       const atom * Otheratom = (*bonditer)->GetOtherAtom(*atomiter);
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/atom.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -72,9 +72,4 @@
 {
   removeFromMolecule();
-  for(BondList::iterator iter=ListOfBonds.begin(); iter!=ListOfBonds.end();){
-    // deleting the bond will invalidate the iterator !!!
-    bond *bond =*(iter++);
-    delete(bond);
-  }
 };
 
@@ -294,4 +289,5 @@
   if (ComponentNr != NULL)
     delete[](ComponentNr);
+  const BondList& ListOfBonds = getListOfBonds();
   ComponentNr = new int[ListOfBonds.size()+1];
   for (int i=ListOfBonds.size()+1;i--;)
Index: src/atom_bondedparticle.cpp
===================================================================
--- src/atom_bondedparticle.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/atom_bondedparticle.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -27,4 +27,5 @@
 #include "CodePatterns/Log.hpp"
 #include "CodePatterns/Verbose.hpp"
+#include "WorldTime.hpp"
 
 /** Constructor of class BondedParticle.
@@ -32,4 +33,5 @@
 BondedParticle::BondedParticle()
 {
+  ListOfBonds.push_back(BondList());
 };
 
@@ -38,8 +40,13 @@
 BondedParticle::~BondedParticle()
 {
-  BondList::const_iterator Runner;
-  while (!ListOfBonds.empty()) {
-    Runner = ListOfBonds.begin();
-    removewithoutcheck(*Runner);
+  BondList::iterator Runner;
+  for (std::vector<BondList>::iterator iter = ListOfBonds.begin();
+      !ListOfBonds.empty();
+      iter = ListOfBonds.begin()) {
+    while (!(*iter).empty()) {
+      Runner = (*iter).begin();
+      removewithoutcheck(*Runner);
+    }
+    ListOfBonds.erase(iter);
   }
 };
@@ -58,4 +65,5 @@
 void BondedParticle::OutputBondOfAtom() const
 {
+  const BondList& ListOfBonds = getListOfBonds();
   DoLog(4) && (Log() << Verbose(4) << "Atom " << getName() << "/" << nr << " with " << ListOfBonds.size() << " bonds: " << endl);
   int TotalDegree = 0;
@@ -72,4 +80,5 @@
 void BondedParticle::OutputAdjacency(ofstream * const AdjacencyFile) const
 {
+  const BondList& ListOfBonds = getListOfBonds();
   *AdjacencyFile << nr << "\t";
   for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
@@ -84,4 +93,5 @@
 void BondedParticle::OutputBonds(ofstream * const BondFile) const
 {
+  const BondList& ListOfBonds = getListOfBonds();
   for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
     if (nr < (*Runner)->GetOtherAtom(this)->nr)
@@ -113,4 +123,5 @@
   if (Binder != NULL) {
     if (Binder->Contains(this)) {
+      BondList& ListOfBonds = getListOfBonds();
       ListOfBonds.push_back(Binder);
       status = true;
@@ -125,4 +136,5 @@
 
 /** Removes a given bond from atom::ListOfBonds.
+ * @param _step time step to access
  * \param *Binder bond to remove
  */
@@ -132,4 +144,5 @@
   if (Binder != NULL) {
     if (Binder->Contains(this)) {
+      BondList& ListOfBonds = getListOfBonds();
       ListOfBonds.remove(Binder);
       status = true;
@@ -148,4 +161,5 @@
 void BondedParticle::UnregisterAllBond()
 {
+  BondList& ListOfBonds = getListOfBonds();
   ListOfBonds.clear();
 };
@@ -165,4 +179,5 @@
   //Log() << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
   if ((int)(getType()->getNoValenceOrbitals()) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
+    const BondList& ListOfBonds = getListOfBonds();
     for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
       OtherWalker = (*Runner)->GetOtherAtom(this);
@@ -170,5 +185,6 @@
       //Log() << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << OtherNoBonds << "?" << endl;
       if ((int)(OtherWalker->getType()->getNoValenceOrbitals()) > OtherNoBonds) { // check if possible candidate
-        if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
+        const BondList& OtherListOfBonds = OtherWalker->getListOfBonds();
+        if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherListOfBonds.size())) { // pick the one with fewer number of bonds first
           CandidateBond = (*Runner);
           //Log() << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
@@ -194,4 +210,5 @@
 {
   int NoBonds = 0;
+  const BondList& ListOfBonds = getListOfBonds();
   for (BondList::const_iterator Runner = ListOfBonds.begin();
       Runner != ListOfBonds.end();
@@ -209,5 +226,6 @@
   bool status = false;
 
-  for (BondList::iterator runner = ListOfBonds.begin(); runner != ListOfBonds.end(); runner++) {
+  const BondList& ListOfBonds = getListOfBonds();
+  for (BondList::const_iterator runner = ListOfBonds.begin(); runner != ListOfBonds.end(); runner++) {
     status = status || ((*runner)->Contains(BondPartner));
   }
Index: src/atom_bondedparticleinfo.cpp
===================================================================
--- src/atom_bondedparticleinfo.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/atom_bondedparticleinfo.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -20,4 +20,8 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include "CodePatterns/Assert.hpp"
+
+#include "WorldTime.hpp"
+
 #include "atom_bondedparticleinfo.hpp"
 
@@ -27,10 +31,52 @@
   AdaptiveOrder(0),
   MaxOrder(false)
-{};
+{}
 
 /** Destructor of class BondedParticleInfo.
  */
 BondedParticleInfo::~BondedParticleInfo()
+{}
+
+const BondList& BondedParticleInfo::getListOfBonds() const
 {
-};
+  ASSERT(WorldTime::getTime() < ListOfBonds.size(),
+      "BondedParticleInfo::getBondsAtStep() - Access out of range: "
+      +toString(WorldTime::getTime())
+      +" not in [0,"+toString(ListOfBonds.size())+").");
+  return ListOfBonds[WorldTime::getTime()];
+}
 
+BondList& BondedParticleInfo::getListOfBonds()
+{
+  const unsigned int size = ListOfBonds.size();
+  ASSERT(WorldTime::getTime() <= size,
+      "BondedParticleInfo::getBondsAtStep() - Access out of range: "
+      +toString(WorldTime::getTime())
+      +" not in [0,"+toString(size)+"].");
+  if (WorldTime::getTime() == size) {
+    ListOfBonds.push_back(BondList());
+  }
+  return ListOfBonds[WorldTime::getTime()];
+}
+
+const BondList& BondedParticleInfo::getListOfBondsAtStep(unsigned int _step) const
+{
+  ASSERT(_step < ListOfBonds.size(),
+      "BondedParticleInfo::getBondsAtStep() - Access out of range: "
+      +toString(WorldTime::getTime())
+      +" not in [0,"+toString(ListOfBonds.size())+").");
+  return ListOfBonds[_step];
+}
+
+BondList& BondedParticleInfo::getListOfBondsAtStep(unsigned int _step)
+{
+  const unsigned int size = ListOfBonds.size();
+  ASSERT(_step <= size,
+      "BondedParticleInfo::getBondsAtStep() - Access out of range: "
+      +toString(_step)
+      +" not in [0,"+toString(size)+"].");
+  if (_step == size) {
+    ListOfBonds.push_back(BondList());
+  }
+  return ListOfBonds[_step];
+}
Index: src/atom_bondedparticleinfo.hpp
===================================================================
--- src/atom_bondedparticleinfo.hpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/atom_bondedparticleinfo.hpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -24,4 +24,5 @@
 
 class bond;
+class BondedParticle;
 
 #define BondList list<bond *>
@@ -30,13 +31,44 @@
 
 class BondedParticleInfo {
+  friend class BondedParticle;
 public:
   unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
   bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
-  BondList ListOfBonds; //!< list of all bonds
 
   BondedParticleInfo();
   virtual ~BondedParticleInfo();
 
+  /** Const accessor to ListOfBonds of WorldTime::CurrentTime.
+   *
+   * @return ListOfBonds[WorldTime::CurrentTime]
+   */
+  const BondList& getListOfBonds() const;
+  /** Accessor to ListOfBonds of WorldTime::CurrentTime.
+   *
+   * Note, new empty BondList is returned if array entry at upper boundary is
+   * accessed. Beyond std will issue exception due to out-of-range access.
+   *
+   * @return ListOfBonds[WorldTime::CurrentTime]
+   */
+  BondList& getListOfBonds();
+
+  /** Const Accessor ListOfBonds of any present time step.
+   *
+   * @param _step time step to access
+   * @return ListOfBonds[_step].
+   */
+  const BondList& getListOfBondsAtStep(unsigned int _step) const;
+  /** Accessor ListOfBonds of any present time step.
+   *
+   * Note, new empty BondList is returned if array entry at upper boundary is
+   * accessed. Beyond std will issue exception due to out-of-range access.
+   *
+   * @param _step time step to access
+   * @return ListOfBonds[_step].
+   */
+  BondList& getListOfBondsAtStep(unsigned int _step);
+
 private:
+  std::vector<BondList> ListOfBonds; //!< list of all bonds
 
 };
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/boundary.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -818,8 +818,13 @@
   const int FillerCount = filler->getAtomCount();
   DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
-  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
-    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner) {
       if ((*BondRunner)->leftatom == *AtomRunner)
         DoLog(2) && (Log() << Verbose(2) << "  " << *(*BondRunner) << endl);
+    }
+  }
 
   atom * CopyAtoms[FillerCount];
@@ -912,6 +917,9 @@
         }
         // go through all bonds and add as well
-        for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
-          for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+        for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
+          const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+          for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+              BondRunner != ListOfBonds.end();
+              ++BondRunner)
             if ((*BondRunner)->leftatom == *AtomRunner) {
               Binder = (*BondRunner);
@@ -921,4 +929,5 @@
               }
             }
+        }
       }
   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
@@ -1084,8 +1093,12 @@
   //const int FillerCount = filler->getAtomCount();
   DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
-  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
-    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner)
         DoLog(2) && (Log() << Verbose(2) << "  " << *(*BondRunner) << endl);
+  }
 
   // calculate filler grid in [0,1]^3
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/config.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -965,5 +965,5 @@
   int MaxNeighbours = 0;
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
-    const int count = (*iter)->ListOfBonds.size();
+    const int count = (*iter)->getListOfBonds().size();
     if (MaxNeighbours < count)
       MaxNeighbours = count;
@@ -979,7 +979,11 @@
     *output << static_cast<double>((*iter)->getType()->getValence()) << "\t";
     *output << (*iter)->getType()->getSymbol() << "\t";
-    for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
+    const BondList& ListOfBonds = (*iter)->getListOfBonds();
+    for (BondList::const_iterator runner = ListOfBonds.begin();
+        runner != ListOfBonds.end();
+        runner++) {
       *output << (*runner)->GetOtherAtom(*iter)->nr << "\t";
-    for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
+    }
+    for(int i = ListOfBonds.size(); i < MaxNeighbours; i++)
       *output << "-\t";
     *output << endl;
@@ -1017,5 +1021,5 @@
   for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
-      const int count = (*iter)->ListOfBonds.size();
+      const int count = (*iter)->getListOfBonds().size();
       if (MaxNeighbours < count)
         MaxNeighbours = count;
@@ -1051,7 +1055,11 @@
         *output << (double)(*iter)->getType()->getValence() << "\t";
         *output << (*iter)->getType()->getSymbol() << "\t";
-        for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
+        const BondList& ListOfBonds = (*iter)->getListOfBonds();
+        for (BondList::const_iterator runner = ListOfBonds.begin();
+            runner != ListOfBonds.end();
+            runner++) {
           *output << LocalNotoGlobalNoMap[ (*runner)->GetOtherAtom((*iter))->getId() ] << "\t";
-        for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
+        }
+        for(int i = ListOfBonds.size(); i < MaxNeighbours; i++)
           *output << "-\t";
         *output << endl;
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/molecule.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -26,27 +26,26 @@
 #include <gsl/gsl_heapsort.h>
 
-#include "World.hpp"
 #include "atom.hpp"
 #include "bond.hpp"
+#include "Box.hpp"
+#include "CodePatterns/enumeration.hpp"
+#include "CodePatterns/Log.hpp"
 #include "config.hpp"
 #include "element.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 #include "graph.hpp"
 #include "Helpers/helpers.hpp"
 #include "LinearAlgebra/leastsquaremin.hpp"
+#include "LinearAlgebra/Plane.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
 #include "linkedcell.hpp"
 #include "lists.hpp"
-#include "CodePatterns/Log.hpp"
 #include "molecule.hpp"
-
 #include "periodentafel.hpp"
 #include "tesselation.hpp"
-#include "LinearAlgebra/Vector.hpp"
-#include "LinearAlgebra/RealSpaceMatrix.hpp"
 #include "World.hpp"
-#include "Box.hpp"
-#include "LinearAlgebra/Plane.hpp"
-#include "Exceptions/LinearDependenceException.hpp"
-
-#include "CodePatterns/enumeration.hpp"
+#include "WorldTime.hpp"
+
 
 /************************************* Functions for class molecule *********************************/
@@ -383,15 +382,20 @@
       break;
     case 2:
-      // determine two other bonds (warning if there are more than two other) plus valence sanity check
-      for (BondList::const_iterator Runner = TopOrigin->ListOfBonds.begin(); Runner != TopOrigin->ListOfBonds.end(); (++Runner)) {
-        if ((*Runner) != TopBond) {
-          if (FirstBond == NULL) {
-            FirstBond = (*Runner);
-            FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
-          } else if (SecondBond == NULL) {
-            SecondBond = (*Runner);
-            SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
-          } else {
-            DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->getName());
+      {
+        // determine two other bonds (warning if there are more than two other) plus valence sanity check
+        const BondList& ListOfBonds = TopOrigin->getListOfBonds();
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != ListOfBonds.end();
+            ++Runner) {
+          if ((*Runner) != TopBond) {
+            if (FirstBond == NULL) {
+              FirstBond = (*Runner);
+              FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
+            } else if (SecondBond == NULL) {
+              SecondBond = (*Runner);
+              SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
+            } else {
+              DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->getName());
+            }
           }
         }
@@ -642,6 +646,9 @@
 
   // copy all bonds
-  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner) {
         bond *Binder = (*BondRunner);
@@ -660,4 +667,5 @@
         NewBond->Type = Binder->Type;
       }
+  }
   // correct fathers
   for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
@@ -749,6 +757,7 @@
   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
   BondList::const_iterator ForeRunner;
-  while (!BondPartner->ListOfBonds.empty()) {
-    ForeRunner = BondPartner->ListOfBonds.begin();
+  BondList& ListOfBonds = BondPartner->getListOfBonds();
+  while (!ListOfBonds.empty()) {
+    ForeRunner = ListOfBonds.begin();
     RemoveBond(*ForeRunner);
   }
Index: src/molecule_fragmentation.cpp
===================================================================
--- src/molecule_fragmentation.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/molecule_fragmentation.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -56,5 +56,6 @@
   // get maximum bond degree
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c;
+    const BondList& ListOfBonds = (*iter)->getListOfBonds();
+    c = (ListOfBonds.size() > c) ? ListOfBonds.size() : c;
   }
   FragmentCount = NoNonHydrogen*(1 << (c*order));
@@ -670,6 +671,7 @@
         atomiter != (*iter)->end();
         ++atomiter) {
-      for (BondList::iterator bonditer = (*atomiter)->ListOfBonds.begin();
-          bonditer != (*atomiter)->ListOfBonds.end();
+      const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
+      for (BondList::const_iterator bonditer = ListOfBonds.begin();
+          bonditer != ListOfBonds.end();
           ++bonditer) {
         (*bonditer)->BondDegree = 1;
@@ -1009,5 +1011,8 @@
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
-      for (BondList::const_iterator BondRunner = FatherOfRunner->ListOfBonds.begin(); BondRunner != FatherOfRunner->ListOfBonds.end(); (++BondRunner)) {
+      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
+      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+          BondRunner != ListOfBonds.end();
+          ++BondRunner) {
         OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
 //        Log() << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
@@ -1416,5 +1421,8 @@
       // go through all its bonds
       DoLog(1) && (Log() << Verbose(1) << "Going through all bonds of Walker." << endl);
-      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      const BondList& ListOfBonds = Walker->getListOfBonds();
+      for (BondList::const_iterator Runner = ListOfBonds.begin();
+          Runner != ListOfBonds.end();
+          ++Runner) {
         OtherWalker = (*Runner)->GetOtherAtom(Walker);
         if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
@@ -1811,6 +1819,9 @@
     // scan all bonds
     flag = false;
-    for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner)
-      for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); (!flag) && (BondRunner != (*AtomRunner)->ListOfBonds.end()); ++BondRunner) {
+    for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner) {
+      const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+      for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+          (!flag) && (BondRunner != ListOfBonds.end());
+          ++BondRunner) {
         Binder = (*BondRunner);
         for (int i=NDIM;i--;) {
@@ -1832,4 +1843,5 @@
         }
       }
+    }
     //if (flag) {
     if (0) {
@@ -1853,5 +1865,8 @@
         ColorList[Walker->nr] = black;    // mark as explored
         *Walker += Translationvector; // translate
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        const BondList& ListOfBonds = Walker->getListOfBonds();
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != ListOfBonds.end();
+            ++Runner) {
           if ((*Runner) != Binder) {
             OtherWalker = (*Runner)->GetOtherAtom(Walker);
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/molecule_geometry.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -384,5 +384,8 @@
         Testvector = inversematrix * (*iter)->getPosition();
         Translationvector.Zero();
-        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
+        const BondList& ListOfBonds = (*iter)->getListOfBonds();
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != 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++) {
@@ -404,5 +407,7 @@
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
-        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != ListOfBonds.end();
+            ++Runner) {
           if ((*Runner)->GetOtherAtom((*iter))->getType()->getAtomicNumber() == 1) {
             Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->getPosition();
Index: src/molecule_graph.cpp
===================================================================
--- src/molecule_graph.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/molecule_graph.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -35,4 +35,5 @@
 #include "molecule.hpp"
 #include "World.hpp"
+#include "WorldTime.hpp"
 #include "Box.hpp"
 
@@ -149,8 +150,12 @@
   DoLog(0) && (Log() << Verbose(0) << "Begin of CreateAdjacencyList." << endl);
   // remove every bond from the list
-  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
+  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::iterator BondRunner = ListOfBonds.begin();
+        !ListOfBonds.empty();
+        BondRunner = ListOfBonds.begin())
       if ((*BondRunner)->leftatom == *AtomRunner)
         delete((*BondRunner));
+  }
   BondCount = 0;
 
@@ -239,7 +244,9 @@
 bool molecule::hasBondStructure() const
 {
-  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    if (!(*AtomRunner)->ListOfBonds.empty())
+  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    if (!ListOfBonds.empty())
       return true;
+  }
   return false;
 }
@@ -251,8 +258,12 @@
 {
   unsigned int counter = 0;
-  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner)
         counter++;
+  }
   return counter;
 }
@@ -264,9 +275,13 @@
 {
   DoLog(1) && (Log() << Verbose(1) << endl << "From contents of bond chain list:");
-  for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner)
-    for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner) {
         DoLog(0) && (Log() << Verbose(0) << *(*BondRunner) << "\t" << endl);
       }
+  }
   DoLog(0) && (Log() << Verbose(0) << endl);
 }
@@ -315,6 +330,7 @@
   MoleculeLeafClass *Subgraphs = NULL;
   std::deque<bond *> *BackEdgeStack = NULL;
-  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    if ((!(*AtomRunner)->ListOfBonds.empty()) && ((*(*AtomRunner)->ListOfBonds.begin())->Type == Undetermined)) {
+  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    if ((!ListOfBonds.empty()) && ((*ListOfBonds.begin())->Type == Undetermined)) {
       DoLog(0) && (Log() << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl);
       Subgraphs = DepthFirstSearchAnalysis(BackEdgeStack);
@@ -327,9 +343,14 @@
       break;
     }
-  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  }
+  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner)
         if ((*BondRunner)->Cyclic)
           NoCyclicBonds++;
+  }
   delete (BackEdgeStack);
   return NoCyclicBonds;
@@ -637,6 +658,9 @@
 {
   NoCyclicBonds = 0;
-  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner)
         if ((*BondRunner)->rightatom->LowpointNr == (*BondRunner)->leftatom->LowpointNr) { // cyclic ??
@@ -644,4 +668,5 @@
           NoCyclicBonds++;
         }
+  }
 }
 ;
@@ -662,10 +687,12 @@
 void molecule::OutputGraphInfoPerBond() const
 {
-  bond *Binder = NULL;
   DoLog(1) && (Log() << Verbose(1) << "Final graph info for each bond is:" << endl);
-  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner) {
-        Binder = *BondRunner;
+        const bond *Binder = *BondRunner;
         if (DoLog(2)) {
           ostream &out = (Log() << Verbose(2));
@@ -681,4 +708,5 @@
           DoLog(3) && (Log() << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl);
       }
+  }
 }
 ;
@@ -762,5 +790,8 @@
     BFS.BFSStack->pop_front();
     DoLog(2) && (Log() << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *BFS.Root << "." << endl);
-    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+    const BondList& ListOfBonds = Walker->getListOfBonds();
+    for (BondList::const_iterator Runner = ListOfBonds.begin();
+        Runner != ListOfBonds.end();
+        ++Runner) {
       if ((*Runner) != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
         OtherAtom = (*Runner)->GetOtherAtom(Walker);
@@ -895,7 +926,10 @@
     BFS.BFSStack->pop_front();
     //Log() << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+    const BondList& ListOfBonds = Walker->getListOfBonds();
+    for (BondList::const_iterator Runner = ListOfBonds.begin();
+        Runner != ListOfBonds.end();
+        ++Runner) {
       // "removed (*Runner) != BackEdge) || " from next if, is u
-      if ((Walker->ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
+      if ((ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
         OtherAtom = (*Runner)->GetOtherAtom(Walker);
         //Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
@@ -1017,5 +1051,6 @@
   size_t i = 0;
   if (vertex != NULL) {
-    for (; i < vertex->ListOfBonds.size(); i++) {
+    const BondList& ListOfBonds = vertex->getListOfBonds();
+    for (; i < ListOfBonds.size(); i++) {
       if (vertex->ComponentNr[i] == -1) { // check if not yet used
         vertex->ComponentNr[i] = nr;
@@ -1024,5 +1059,5 @@
         break; // breaking here will not cause error!
     }
-    if (i == vertex->ListOfBonds.size()) {
+    if (i == ListOfBonds.size()) {
       DoeLog(0) && (eLog()<< Verbose(0) << "Error: All Component entries are already occupied!" << endl);
       performCriticalExit();
@@ -1041,5 +1076,8 @@
 bond * molecule::FindNextUnused(atom *vertex) const
 {
-  for (BondList::const_iterator Runner = vertex->ListOfBonds.begin(); Runner != vertex->ListOfBonds.end(); (++Runner))
+  const BondList& ListOfBonds = vertex->getListOfBonds();
+  for (BondList::const_iterator Runner = ListOfBonds.begin();
+      Runner != ListOfBonds.end();
+      ++Runner)
     if ((*Runner)->IsUsed() == white)
       return ((*Runner));
@@ -1053,8 +1091,12 @@
 void molecule::ResetAllBondsToUnused() const
 {
-  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
-    for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner)
       if ((*BondRunner)->leftatom == *AtomRunner)
         (*BondRunner)->ResetUsed();
+  }
 }
 ;
@@ -1168,6 +1210,9 @@
 
   //Log() << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
-  if (CurrentBondsOfAtom == Walker->ListOfBonds.size()) {
-    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+  const BondList& ListOfBonds = Walker->getListOfBonds();
+  if (CurrentBondsOfAtom == ListOfBonds.size()) {
+    for (BondList::const_iterator Runner = ListOfBonds.begin();
+        Runner != ListOfBonds.end();
+        ++Runner) {
       id = (*Runner)->GetOtherAtom(Walker)->nr;
       j = 0;
@@ -1185,5 +1230,5 @@
     //Log() << Verbose(0) << endl;
   } else {
-    DoLog(0) && (Log() << Verbose(0) << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << Walker->ListOfBonds.size() << "." << endl);
+    DoLog(0) && (Log() << Verbose(0) << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << ListOfBonds.size() << "." << endl);
     status = false;
   }
@@ -1271,6 +1316,9 @@
   do { // go through all bonds and push local ones
     Walker = ListOfLocalAtoms[Binder->leftatom->nr]; // get one atom in the reference molecule
-    if (Walker != NULL) // if this Walker exists in the subgraph ...
-      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+    if (Walker != NULL) { // if this Walker exists in the subgraph ...
+      const BondList& ListOfBonds = Walker->getListOfBonds();
+      for (BondList::const_iterator Runner = ListOfBonds.begin();
+          Runner != ListOfBonds.end();
+          ++Runner) {
         OtherAtom = (*Runner)->GetOtherAtom(Walker);
         if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->nr]) { // found the bond
@@ -1280,4 +1328,5 @@
         }
       }
+    }
     ASSERT(!ReferenceStack->empty(), "molecule::PickLocalBackEdges() - ReferenceStack is empty!");
     Binder = ReferenceStack->front(); // loop the stack for next item
@@ -1423,6 +1472,9 @@
     Walker = BFS.BFSStack->front(); // pop oldest added
     BFS.BFSStack->pop_front();
-    DoLog(1) && (Log() << Verbose(1) << "Current Walker is: " << Walker->getName() << ", and has " << Walker->ListOfBonds.size() << " bonds." << endl);
-    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+    const BondList& ListOfBonds = Walker->getListOfBonds();
+    DoLog(1) && (Log() << Verbose(1) << "Current Walker is: " << Walker->getName() << ", and has " << ListOfBonds.size() << " bonds." << endl);
+    for (BondList::const_iterator Runner = ListOfBonds.begin();
+        Runner != ListOfBonds.end();
+        ++Runner) {
       if ((*Runner) != NULL) { // don't look at bond equal NULL
         Binder = (*Runner);
@@ -1496,5 +1548,8 @@
         status = false;
       } else {
-        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
+        const BondList& ListOfBonds = (*iter)->getListOfBonds();
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != 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
@@ -1534,5 +1589,4 @@
 
 /** For a given keyset \a *Fragment, checks whether it is connected in the current molecule.
- * \param *out output stream for debugging
  * \param *Fragment Keyset of fragment's vertices
  * \return true - connected, false - disconnected
@@ -1558,5 +1612,8 @@
       for (KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
         Walker2 = FindAtom(*runners);
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        const BondList& ListOfBonds = Walker->getListOfBonds();
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != ListOfBonds.end();
+            ++Runner) {
           if ((*Runner)->GetOtherAtom(Walker) == Walker2) {
             BondStatus = true;
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/moleculelist.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -363,11 +363,11 @@
     // 2. take every hydrogen that is a saturated one
     for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
-      //Log() << Verbose(1) << "(*iter): " << *(*iter) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
+      //Log() << Verbose(1) << "(*iter): " << *(*iter) << " with first bond " << *((*iter)->getListOfBonds().begin()) << "." << endl;
       if (((*iter)->getType()->getAtomicNumber() == 1) && (((*iter)->father == NULL)
           || ((*iter)->father->getType()->getAtomicNumber() != 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;
+          //Log() << Verbose(2) << "Runner: " << *(*runner) << " with first bond " << *((*iter)->getListOfBonds().begin()) << "." << endl;
           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
-          Binder = *((*runner)->ListOfBonds.begin());
+          Binder = *((*runner)->getListOfBonds().begin());
           if (((*runner)->getType()->getAtomicNumber() == 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
@@ -805,13 +805,21 @@
     DoLog(1) && (Log() << Verbose(1) << "Creating adjacency list for subgraph " << Leaf << "." << endl);
     // remove every bond from the list
-    for(molecule::iterator AtomRunner = Leaf->begin(); AtomRunner != Leaf->end(); ++AtomRunner)
-      for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
+    for(molecule::iterator AtomRunner = Leaf->begin(); AtomRunner != Leaf->end(); ++AtomRunner) {
+      BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+      for(BondList::iterator BondRunner = ListOfBonds.begin();
+          !ListOfBonds.empty();
+          BondRunner = ListOfBonds.begin())
         if ((*BondRunner)->leftatom == *AtomRunner)
           delete((*BondRunner));
+    }
+
 
     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)) {
+      const BondList& ListOfBonds = Father->getListOfBonds();
+      for (BondList::const_iterator Runner = ListOfBonds.begin();
+          Runner != ListOfBonds.end();
+          ++Runner) {
         OtherWalker = ListOfLocalAtoms[(*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr]; // local copy of current bond partner of walker
         if (OtherWalker != NULL) {
Index: src/unittests/ListOfBondsUnitTest.cpp
===================================================================
--- src/unittests/ListOfBondsUnitTest.cpp	(revision c0d9ebc6c65d3d9cc5fc8234ff04f147c042c124)
+++ src/unittests/ListOfBondsUnitTest.cpp	(revision 9d83b69ef70cf27f82b8fcaaf764866bcda9a143)
@@ -126,7 +126,10 @@
 
   // check that bond is present in both atoms
-  bond *TestBond1 = *(atom1->ListOfBonds.begin());
+  BondList::const_iterator bonditer;
+  bonditer = atom1->getListOfBonds().begin();
+  bond *TestBond1 = *bonditer;
   CPPUNIT_ASSERT_EQUAL( TestBond1, Binder );
-  bond *TestBond2 = *(atom2->ListOfBonds.begin());
+  bonditer = atom2->getListOfBonds().begin();
+  bond *TestBond2 = *bonditer;
   CPPUNIT_ASSERT_EQUAL( TestBond2, Binder );
 };
@@ -153,6 +156,12 @@
 
   // check if removed from atoms
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom2->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom2->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
 
   // check if removed from molecule
@@ -185,7 +194,16 @@
 
   // check that all are present
-  CPPUNIT_ASSERT_EQUAL( (size_t) 2, atom1->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 2, atom2->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 2, atom3->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 2, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom2->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 2, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom3->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 2, ListOfBonds.size() );
+  }
 
   // remove bond
@@ -193,7 +211,16 @@
 
   // check if removed from atoms
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 1, atom2->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 1, atom3->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom2->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 1, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom3->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 1, ListOfBonds.size() );
+  }
 
   // check if removed from molecule
@@ -223,6 +250,12 @@
 
   // check if removed from atoms
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom2->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom2->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
 
   // check if removed from molecule
@@ -251,5 +284,8 @@
 
   // check bond if removed from other atom
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
 
   // check if removed from molecule
@@ -278,6 +314,12 @@
   CPPUNIT_ASSERT( Binder != NULL );
 
-  CPPUNIT_ASSERT_EQUAL( (size_t) 1, atom1->ListOfBonds.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t) 1, atom2->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 1, ListOfBonds.size() );
+  }
+  {
+    const BondList& ListOfBonds = atom2->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 1, ListOfBonds.size() );
+  }
 
   CPPUNIT_ASSERT_EQUAL( true, TestMolecule->hasBondStructure() );
@@ -287,5 +329,8 @@
 
   // check bond if removed from other atom
-  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
+  {
+    const BondList& ListOfBonds = atom1->getListOfBonds();
+    CPPUNIT_ASSERT_EQUAL( (size_t) 0, ListOfBonds.size() );
+  }
 
   // check if removed from molecule
