Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/atom.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -9,4 +9,5 @@
 #include "config.hpp"
 #include "element.hpp"
+#include "lists.hpp"
 #include "memoryallocator.hpp"
 #include "parser.hpp"
@@ -67,4 +68,10 @@
 atom::~atom()
 {
+  BondList::const_iterator Runner;
+  while (!ListOfBonds.empty()) {
+    Runner = ListOfBonds.begin();
+    removewithoutcheck(*Runner);
+  }
+  unlink(this);
   Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
   Free<char>(&Name, "atom::~atom: *Name");
@@ -118,4 +125,15 @@
 {
   return (node->IsInParallelepiped(offset, parallelepiped));
+};
+
+/** Counts the number of bonds weighted by bond::BondDegree.
+ * \param bonds times bond::BondDegree
+ */
+int atom::CountBonds() const
+{
+  int NoBonds = 0;
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    NoBonds += (*Runner)->BondDegree;
+  return NoBonds;
 };
 
@@ -215,11 +233,32 @@
 };
 
-/** Prints all bonds of this atom from given global lists.
+/** Output graph info of this atom.
+ * \param *out output stream
+ */
+void atom::OutputGraphInfo(ofstream *out) const
+{
+  *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
+  OutputComponentNumber(out);
+  *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
+};
+
+/** Output a list of flags, stating whether the bond was visited or not.
+ * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
+ * \param *out output stream for debugging
+ */
+void atom::OutputComponentNumber(ofstream *out) const
+{
+  if (ComponentNr != NULL) {
+    for (int i=0; ComponentNr[i] != -1; i++)
+      *out << ComponentNr[i] << " ";
+  }
+};
+
+
+/** Prints all bonds of this atom with total degree.
  * \param *out stream to output to
- * \param *NumberOfBondsPerAtom array with number of bonds per atomic index
- * \param ***ListOfBondsPerAtom array per atomic index of array with pointer to bond
  * \return true - \a *out present, false - \a *out is NULL
  */
-bool atom::OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const
+bool atom::OutputBondOfAtom(ofstream *out) const
 {
   if (out != NULL) {
@@ -227,9 +266,9 @@
     if (type->Z != 1) {   // regard only non-hydrogen
 #endif
-      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: ";
+      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
       int TotalDegree = 0;
-      for (int j=0;j<NumberOfBondsPerAtom[nr];j++) {
-        *out << *ListOfBondsPerAtom[nr][j] << "\t";
-        TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree;
+      for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
+        *out << **Runner << "\t";
+        TotalDegree += (*Runner)->BondDegree;
       }
       *out << " -- TotalDegree: " << TotalDegree << endl;
@@ -240,4 +279,97 @@
   } else
     return false;
+};
+
+/** Output of atom::nr along with all bond partners.
+ * \param *AdjacencyFile output stream
+ */
+void atom::OutputAdjacency(ofstream *AdjacencyFile) const
+{
+  *AdjacencyFile << nr << "\t";
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
+  *AdjacencyFile << endl;
+};
+
+/** Puts a given bond into atom::ListOfBonds.
+ * \param *Binder bond to insert
+ */
+bool atom::RegisterBond(bond *Binder)
+{
+  bool status = false;
+  if (Binder != NULL) {
+    if (Binder->Contains(this)) {
+      ListOfBonds.push_back(Binder);
+      status = true;
+    } else {
+      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
+    }
+  } else {
+    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
+  }
+  return status;
+};
+
+/** Removes a given bond from atom::ListOfBonds.
+ * \param *Binder bond to remove
+ */
+bool atom::UnregisterBond(bond *Binder)
+{
+  bool status = false;
+  if (Binder != NULL) {
+    if (Binder->Contains(this)) {
+      ListOfBonds.remove(Binder);
+      status = true;
+    } else {
+      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
+    }
+  } else {
+    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
+  }
+  return status;
+};
+
+/** Removes all bonds from atom::ListOfBonds.
+ * \note Does not do any memory de-allocation.
+ */
+void atom::UnregisterAllBond()
+{
+  ListOfBonds.clear();
+};
+
+/** Corrects the bond degree by one at most if necessary.
+ * \param *out output stream for debugging
+ */
+int atom::CorrectBondDegree(ofstream *out)
+{
+  int NoBonds = 0;
+  int OtherNoBonds = 0;
+  int FalseBondDegree = 0;
+  atom *OtherWalker = NULL;
+  bond *CandidateBond = NULL;
+
+  *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+  NoBonds = CountBonds();
+  if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
+    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
+      OtherWalker = (*Runner)->GetOtherAtom(this);
+      OtherNoBonds = OtherWalker->CountBonds();
+      *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+      if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
+        if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
+          CandidateBond = (*Runner);
+          *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
+        }
+      }
+    }
+    if ((CandidateBond != NULL)) {
+      CandidateBond->BondDegree++;
+      *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
+    } else {
+      *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
+      FalseBondDegree++;
+    }
+  }
+  return FalseBondDegree;
 };
 
@@ -345,5 +477,5 @@
 {
   *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
-  //cout << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << "." << endl;
+  //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
 }
 
Index: src/atom.hpp
===================================================================
--- src/atom.hpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/atom.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -19,4 +19,5 @@
 
 #include <iostream>
+#include <list>
 #include <vector>
 
@@ -32,4 +33,6 @@
 class ForceMatrix;
 class Vector;
+
+#define BondList list<bond *>
 
 /********************************************** declarations *******************************/
@@ -50,4 +53,5 @@
     Vector v;       //!< velocity vector of atom, giving last velocity within cell
     Vector F;       //!< Force vector of atom, giving last force within cell
+    BondList ListOfBonds; //!< list of all bonds
     element *type;  //!< pointing to element
     atom *previous; //!< previous atom in molecule list
@@ -76,5 +80,7 @@
   bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
   bool OutputTrajectoryXYZ(ofstream *out, int step) const;
-  bool OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const;
+  bool OutputBondOfAtom(ofstream *out) const;
+  void OutputAdjacency(ofstream *AdjacencyFile) const;
+
 
   void EqualsFather ( atom *ptr, atom **res );
@@ -96,4 +102,12 @@
   // bond order stuff
   void OutputOrder(ofstream *file);
+  void OutputGraphInfo(ofstream *out) const;
+  void OutputComponentNumber(ofstream *out) const;
+  int CountBonds() const;
+  int CorrectBondDegree(ofstream *out);
+
+  bool RegisterBond(bond *Binder);
+  bool UnregisterBond(bond *Binder);
+  void UnregisterAllBond();
 
   // constraint potential and dynamics stuff
Index: src/bond.cpp
===================================================================
--- src/bond.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/bond.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -53,4 +53,5 @@
   Cyclic = false;
 };
+
 bond::bond(atom *left, atom *right)
 {
@@ -77,10 +78,9 @@
 {
   // remove this node from the list structure
-  if (previous != NULL) {
-    previous->next = next;
-  }
-  if (next != NULL) {
-    next->previous = previous;
-  }
+  if (leftatom != NULL)
+    leftatom->UnregisterBond(this);
+  if (rightatom != NULL)
+  rightatom->UnregisterBond(this);
+  unlink(this);
 };
 
@@ -95,5 +95,5 @@
  * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
  */
-atom * bond::GetOtherAtom(atom *Atom) const
+atom * bond::GetOtherAtom(const atom *Atom) const
 {
   if(leftatom == Atom) 
@@ -105,21 +105,4 @@
 };
 
-/** Get the other atom in a bond if one is specified.
- * \param *Atom the pointer to the one atom
- * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
- */
-bond * bond::GetFirstBond()
-{
-  return GetFirst(this);
-};
-
-/** Get the other atom in a bond if one is specified.
- * \param *Atom the pointer to the one atom
- * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
- */
-bond * bond::GetLastBond()
-{
-  return GetLast(this);
-};
 
 /** Returns whether vertex was used in DFS.
Index: src/bond.hpp
===================================================================
--- src/bond.hpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/bond.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -25,7 +25,9 @@
 
 /** Bonds between atoms.
- * Class incorporates bonds between atoms in a molecule,
- * used to derive tge fragments in many-body bond order
- * calculations.
+ * Class incorporates bonds between atoms in a molecule.
+ * Note that we regard bond always as something in a molecule,
+ * as it is the glue making up the connected subgrapgh and
+ * hence the molecule. Thus, bonds belong globally to the molecule
+ * (and are free'd there) and only locally to the atom classs.
  */
 class bond {
@@ -41,7 +43,5 @@
     enum EdgeType Type;//!< whether this is a tree or back edge
 
-  atom * GetOtherAtom(atom *Atom) const;
-  bond * GetFirstBond();
-  bond * GetLastBond();
+  atom * GetOtherAtom(const atom *Atom) const;
 
   bool MarkUsed(enum Shading color);
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/builder.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -895,5 +895,4 @@
           start = clock();
           mol->CreateAdjacencyList((ofstream *)&cout, bonddistance, configuration->GetIsAngstroem());
-          mol->CreateListOfBondsPerAtom((ofstream *)&cout);
           end = clock();
           cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
@@ -1555,5 +1554,4 @@
                 class StackClass<bond *> *LocalBackEdgeStack = NULL;
                 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
-                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
                 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
                 if (Subgraphs != NULL) {
Index: src/linkedcell.cpp
===================================================================
--- src/linkedcell.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/linkedcell.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -203,4 +203,5 @@
 
 /** Checks whether LinkedCell::n[] plus relative offset is each within [0,N[]].
+ * Note that for this check we don't admonish if out of bounds.
  * \param relative[NDIM] relative offset to current cell
  * \return if all in intervals - true, else -false
@@ -211,6 +212,4 @@
   for(int i=0;i<NDIM;i++)
     status = status && ((n[i]+relative[i] >=0) && (n[i]+relative[i] < N[i]));
-  if (!status)
-  cerr << "ERROR: indices are out of bounds!" << endl;
   return status;
 };
Index: src/lists.hpp
===================================================================
--- src/lists.hpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/lists.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -35,4 +35,6 @@
   if (walker->previous != NULL)
     walker->previous->next = walker->next;
+  walker->next = NULL;
+  walker->previous= NULL;
 };
 
@@ -112,12 +114,10 @@
 {
   X *pointer = start->next;
-  X *walker;
+  X *walker = NULL;
   while (pointer != end) { // go through list
     walker = pointer; // mark current
     pointer = pointer->next; // step onward beforehand
     // remove walker
-    unlink(walker);
-    delete(walker);
-    walker = NULL;
+    removewithoutcheck(walker);
   }
   return true;
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/molecule.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -48,6 +48,4 @@
   NoNonHydrogen = 0;
   NoCyclicBonds = 0;
-  ListOfBondsPerAtom = NULL;
-  NumberOfBondsPerAtom = NULL;
   ElementCount = 0;
   for(int i=MAX_ELEMENTS;i--;)
@@ -66,9 +64,4 @@
 molecule::~molecule()
 {
-  if (ListOfBondsPerAtom != NULL)
-    for(int i=AtomCount;i--;)
-      Free(&ListOfBondsPerAtom[i]);
-  Free(&ListOfBondsPerAtom);
-  Free(&NumberOfBondsPerAtom);
   if (TesselStruct != NULL)
     delete(TesselStruct);
@@ -160,11 +153,9 @@
  * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
  * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
- * \param **BondList list of bonds \a *replacement has (necessary to determine plane for double and triple bonds)
- * \param NumBond  number of bonds in \a **BondList
  * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
  * \return number of atoms added, if < bond::BondDegree then something went wrong
  * \todo double and triple bonds splitting (always use the tetraeder angle!)
  */
-bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem)
+bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
 {
   double bondlength;  // bond length of the bond to be replaced/cut
@@ -177,6 +168,6 @@
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
+  double *matrix;
   bond *Binder = NULL;
-  double *matrix;
 
 //  *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -248,12 +239,12 @@
     case 2:
       // determine two other bonds (warning if there are more than two other) plus valence sanity check
-      for (int i=0;i<NumBond;i++) {
-        if (BondList[i] != TopBond) {
+      for (BondList::const_iterator Runner = TopOrigin->ListOfBonds.begin(); Runner != TopOrigin->ListOfBonds.end(); (++Runner)) {
+        if ((*Runner) != TopBond) {
           if (FirstBond == NULL) {
-            FirstBond = BondList[i];
-            FirstOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
+            FirstBond = (*Runner);
+            FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else if (SecondBond == NULL) {
-            SecondBond = BondList[i];
-            SecondOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
+            SecondBond = (*Runner);
+            SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else {
             *out << Verbose(3) << "WARNING: Detected more than four bonds for atom " << TopOrigin->Name;
@@ -523,5 +514,4 @@
   if (first->next != last) {  // if adjaceny list is present
     copy->BondDistance = BondDistance;
-    copy->CreateListOfBondsPerAtom((ofstream *)&cout);
   }
 
@@ -557,4 +547,6 @@
   if ((atom1 != NULL) && (FindAtom(atom1->nr) != NULL) && (atom2 != NULL) && (FindAtom(atom2->nr) != NULL)) {
     Binder = new bond(atom1, atom2, degree, BondCount++);
+    atom1->RegisterBond(Binder);
+    atom2->RegisterBond(Binder);
     if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
       NoNonBonds++;
@@ -585,5 +577,10 @@
 bool molecule::RemoveBonds(atom *BondPartner)
 {
-  cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  BondList::const_iterator ForeRunner;
+  while (!BondPartner->ListOfBonds.empty()) {
+    ForeRunner = BondPartner->ListOfBonds.begin();
+    RemoveBond(*ForeRunner);
+  }
   return false;
 };
@@ -635,4 +632,5 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
+  RemoveBonds(pointer);
   return remove(pointer, start, end);
 };
@@ -661,5 +659,5 @@
 bool molecule::CleanupMolecule()
 {
-  return (cleanup(start,end) && cleanup(first,last));
+  return (cleanup(first,last) && cleanup(start,end));
 };
 
@@ -773,11 +771,11 @@
 };
 
-/** Outputs contents of molecule::ListOfBondsPerAtom.
+/** Outputs contents of each atom::ListOfBonds.
  * \param *out output stream
  */
 void molecule::OutputListOfBonds(ofstream *out) const
 {
-  *out << Verbose(2) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
-  ActOnAllAtoms (&atom::OutputBondOfAtom, out, NumberOfBondsPerAtom, ListOfBondsPerAtom);
+  *out << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
+  ActOnAllAtoms (&atom::OutputBondOfAtom, out);
   *out << endl;
 };
@@ -900,59 +898,4 @@
   }
   configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
-};
-
-
-/** Creates an 2d array of pointer with an entry for each atom and each bond it has.
- * Updates molecule::ListOfBondsPerAtom, molecule::NumberOfBondsPerAtom by parsing through
- * bond chain list, using molecule::AtomCount and molecule::BondCount.
- * Allocates memory, fills the array and exits
- * \param *out output stream for debugging
- */
-void molecule::CreateListOfBondsPerAtom(ofstream *out)
-{
-  bond *Binder = NULL;
-  *out << Verbose(1) << "Begin of Creating ListOfBondsPerAtom: AtomCount = " << AtomCount << "\tBondCount = " << BondCount << "\tNoNonBonds = " << NoNonBonds << "." << endl;
-
-  // re-allocate memory
-  *out << Verbose(2) << "(Re-)Allocating memory." << endl;
-  if (ListOfBondsPerAtom != NULL) {
-    for(int i=AtomCount;i--;)
-      Free(&ListOfBondsPerAtom[i]);
-    Free(&ListOfBondsPerAtom);
-  }
-  if (NumberOfBondsPerAtom != NULL)
-    Free(&NumberOfBondsPerAtom);
-  ListOfBondsPerAtom = Malloc<bond**>(AtomCount, "molecule::CreateListOfBondsPerAtom: ***ListOfBondsPerAtom");
-  NumberOfBondsPerAtom = Malloc<int>(AtomCount, "molecule::CreateListOfBondsPerAtom: *NumberOfBondsPerAtom");
-
-  // reset bond counts per atom
-  for(int i=AtomCount;i--;)
-    NumberOfBondsPerAtom[i] = 0;
-  // count bonds per atom
-  Binder = first;
-  while (Binder->next != last) {
-    Binder = Binder->next;
-    NumberOfBondsPerAtom[Binder->leftatom->nr]++;
-    NumberOfBondsPerAtom[Binder->rightatom->nr]++;
-  }
-  for(int i=AtomCount;i--;) {
-    // allocate list of bonds per atom
-    ListOfBondsPerAtom[i] = Malloc<bond*>(NumberOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: **ListOfBondsPerAtom[]");
-    // clear the list again, now each NumberOfBondsPerAtom marks current free field
-    NumberOfBondsPerAtom[i] = 0;
-  }
-  // fill the list
-  Binder = first;
-  while (Binder->next != last) {
-    Binder = Binder->next;
-    ListOfBondsPerAtom[Binder->leftatom->nr][NumberOfBondsPerAtom[Binder->leftatom->nr]++] = Binder;
-    ListOfBondsPerAtom[Binder->rightatom->nr][NumberOfBondsPerAtom[Binder->rightatom->nr]++] = Binder;
-  }
-
-  // output list for debugging
-  *out << Verbose(3) << "ListOfBondsPerAtom for each atom:" << endl;
-  ActOnAllAtoms( &atom::OutputBondOfAtom, out, NumberOfBondsPerAtom, ListOfBondsPerAtom );
-
-  *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
 };
 
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/molecule.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -86,7 +86,5 @@
     bond *first;        //!< start of bond list
     bond *last;         //!< end of bond list
-    bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
     int MDSteps;        //!< The number of MD steps in Trajectories
-    int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
     int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
     int BondCount;          //!< number of atoms, brought up-to-date by CountBonds()
@@ -190,4 +188,13 @@
   template <typename T> void SetAtomValueToIndexedArray ( T *array, int element::*index, T atom::*value );
 
+  template <typename res> res SumPerAtom(res (atom::*f)() );
+  template <typename res> res SumPerAtom(res (atom::*f)() const );
+  template <typename res> res SumPerAtom(res (atom::*f)() ) const;
+  template <typename res> res SumPerAtom(res (atom::*f)() const ) const;
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) , T t );
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) const, T t );
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) , T t ) const;
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) const, T t ) const;
+
   /// remove atoms from molecule.
   bool AddAtom(atom *pointer);
@@ -199,5 +206,5 @@
   atom * AddCopyAtom(atom *pointer);
   bool AddXYZFile(string filename);
-  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
+  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem);
   bond * AddBond(atom *first, atom *second, int degree = 1);
   bool RemoveBond(bond *pointer);
@@ -247,8 +254,9 @@
   void CreateAdjacencyListFromDbondFile(ofstream *out, ifstream *output);
   void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
-  void CreateListOfBondsPerAtom(ofstream *out);
   int CorrectBondDegree(ofstream *out);
   void OutputBondsList(ofstream *out);
-  int CountAtomsBonds(int nr);
+  void CyclicBondAnalysis();
+  void OutputGraphInfoPerAtom(ofstream *out);
+  void OutputGraphInfoPerBond(ofstream *out);
 
 
@@ -260,5 +268,4 @@
   void SetNextComponentNumber(atom *vertex, int nr);
   void InitComponentNumbers();
-  void OutputComponentNumber(ofstream *out, atom *vertex);
   void ResetAllBondsToUnused();
   void ResetAllAtomNumbers();
@@ -266,4 +273,6 @@
   bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
   string GetColor(enum Shading color);
+  bond * CopyBond(atom *left, atom *right, bond *CopyBond);
+
 
   molecule *CopyMolecule();
Index: src/molecule_fragmentation.cpp
===================================================================
--- src/molecule_fragmentation.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/molecule_fragmentation.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -31,5 +31,5 @@
 int molecule::GuesstimateFragmentCount(ofstream *out, int order)
 {
-  int c = 0;
+  size_t c = 0;
   int FragmentCount;
   // get maximum bond degree
@@ -37,5 +37,5 @@
   while (Walker->next != end) {
     Walker = Walker->next;
-    c = (NumberOfBondsPerAtom[Walker->nr] > c) ? NumberOfBondsPerAtom[Walker->nr] : c;
+    c = (Walker->ListOfBonds.size() > c) ? Walker->ListOfBonds.size() : c;
   }
   FragmentCount = NoNonHydrogen*(1 << (c*order));
@@ -564,7 +564,4 @@
 
   // ===== 1. Check whether bond structure is same as stored in files ====
-
-  // fill the adjacency list
-  CreateListOfBondsPerAtom(out);
 
   // create lookup table for Atom::nr
@@ -630,5 +627,5 @@
       MolecularWalker = MolecularWalker->next;
       *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
-      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output ListOfBondsPerAtom for debugging
+      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output atom::ListOfBonds for debugging
       if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
         // call BOSSANOVA method
@@ -879,6 +876,6 @@
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
-      for (int i=0;i<mol->NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
-        OtherFather = mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
+      for (BondList::const_iterator BondRunner = FatherOfRunner->ListOfBonds.begin(); BondRunner != FatherOfRunner->ListOfBonds.end(); (++BondRunner)) {
+        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
 //        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
         if (SonList[OtherFather->nr] != NULL) {
@@ -887,5 +884,5 @@
 //            *out << Verbose(3) << "Adding Bond: ";
 //            *out <<
-            Leaf->AddBond(Runner, SonList[OtherFather->nr], mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
+            Leaf->AddBond(Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree);
 //            *out << "." << endl;
             //NumBonds[Runner->nr]++;
@@ -898,8 +895,8 @@
 #ifdef ADDHYDROGEN
           //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          if(!Leaf->AddHydrogenReplacementAtom(out, mol->ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, mol->ListOfBondsPerAtom[FatherOfRunner->nr],mol->NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
+          if(!Leaf->AddHydrogenReplacementAtom(out, (*BondRunner), Runner, FatherOfRunner, OtherFather, IsAngstroem))
             exit(1);
 #endif
-          //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
+          //NumBonds[Runner->nr] += Binder->BondDegree;
         }
       }
@@ -936,5 +933,4 @@
   CreateInducedSubgraphOfFragment(out, this, Leaf, SonList, IsAngstroem);
 
-  Leaf->CreateListOfBondsPerAtom(out);
   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
   Free(&SonList);
@@ -981,5 +977,5 @@
     if (bit) {  // if bit is set, we add this bond partner
       OtherWalker = BondsSet[j]->rightatom;  // rightatom is always the one more distant, i.e. the one to add
-      //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][CurrentCombination] << ", checking on " << *OtherWalker << "." << endl;
+      //*out << Verbose(1+verbosity) << "Current Bond is " << BondsSet[j] << ", checking on " << *OtherWalker << "." << endl;
       *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
       TestKeySetInsert = FragmentSet->insert(OtherWalker->nr);
@@ -1089,7 +1085,7 @@
   NumCombinations = 1 << SetDimension;
 
-  // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
-  // von Endstuecken (aus den Bonds) hinzugefuegt werden und fuer verbleibende ANOVAOrder
-  // rekursiv GraphCrawler in der naechsten Ebene aufgerufen werden
+  // here for all bonds of Walker all combinations of end pieces (from the bonds)
+  // have to be added and for the remaining ANOVA order GraphCrawler be called
+  // recursively for the next level
 
   *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
@@ -1258,6 +1254,6 @@
   atom *OtherWalker = NULL;
   atom *Predecessor = NULL;
+  bond *CurrentEdge = NULL;
   bond *Binder = NULL;
-  bond *CurrentEdge = NULL;
   int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
   int RemainingWalkers = -1;
@@ -1285,7 +1281,6 @@
       // go through all its bonds
       *out << Verbose(1) << "Going through all bonds of Walker." << endl;
-      for (int i=0;i<mol->NumberOfBondsPerAtom[AtomKeyNr];i++) {
-        Binder = mol->ListOfBondsPerAtom[AtomKeyNr][i];
-        OtherWalker = Binder->GetOtherAtom(Walker);
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        OtherWalker = (*Runner)->GetOtherAtom(Walker);
         if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
   #ifdef ADDHYDROGEN
@@ -1293,5 +1288,5 @@
   #endif
                                                               ) {  // skip hydrogens and restrict to fragment
-          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
+          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *(*Runner) << "." << endl;
           // set the label if not set (and push on root stack as well)
           if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
@@ -1712,7 +1707,7 @@
         ColorList[Walker->nr] = black;    // mark as explored
         Walker->x.AddVector(&Translationvector); // translate
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {  // go through all binding partners
-          if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
-            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner) != Binder) {
+            OtherWalker = (*Runner)->GetOtherAtom(Walker);
             if (ColorList[OtherWalker->nr] == white) {
               AtomStack->Push(OtherWalker); // push if yet unexplored
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/molecule_geometry.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -258,5 +258,4 @@
 {
   atom *Walker = start;
-  bond *Binder = NULL;
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   double tmp;
@@ -275,12 +274,11 @@
         Testvector.InverseMatrixMultiplication(matrix);
         Translationvector.Zero();
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          Binder = ListOfBondsPerAtom[Walker->nr][i];
-          if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
+        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 (int j=0;j<NDIM;j++) {
-              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
+              tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
-                cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl;
+                cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl;
                 if (tmp > 0)
                   Translationvector.x[j] -= 1.;
@@ -298,8 +296,7 @@
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          Binder = ListOfBondsPerAtom[Walker->nr][i];
-          if (Binder->GetOtherAtom(Walker)->type->Z == 1) {
-            Testvector.CopyVector(&Binder->GetOtherAtom(Walker)->x);
+        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);
             Testvector.InverseMatrixMultiplication(matrix);
             Testvector.AddVector(&Translationvector);
Index: src/molecule_graph.cpp
===================================================================
--- src/molecule_graph.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/molecule_graph.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -48,6 +48,4 @@
     AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
   }
-
-  CreateListOfBondsPerAtom(out);
 };
 
@@ -64,5 +62,4 @@
  *  -# put each atom into its corresponding cell
  *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
- *  -# create the list of bonds via CreateListOfBondsPerAtom()
  *  -# correct the bond degree iteratively (single->double->triple bond)
  *  -# finally print the bond list to \a *out if desired
@@ -148,12 +145,9 @@
     *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl;
 
-    // create the adjacency list per atom
-    CreateListOfBondsPerAtom(out);
-
     // correct bond degree by comparing valence and bond degree
     CorrectBondDegree(out);
 
     // output bonds for debugging (if bond chain list was correctly installed)
-    OutputBondsList(out);
+    ActOnAllAtoms( &atom::OutputBondOfAtom, out );
   } else
     *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
@@ -173,16 +167,4 @@
   }
   *out << endl;
-};
-
-/** Sums up the number of bonds times bond::BondDegree the atom with index \a nr has.
- * \param nr indexof atom
- * \return number of bonds each weighted with its bond::BondDegree
- */
-int molecule::CountAtomsBonds(int nr)
-{
-  int NoBonds = 0;
-  for(int j=0;j<NumberOfBondsPerAtom[nr];j++)
-    NoBonds += ListOfBondsPerAtom[nr][j]->BondDegree;
-  return NoBonds;
 };
 
@@ -198,46 +180,9 @@
 {
   int No = 0;
-  int NoBonds = 0;
-  int CandidateBondNo = 0;
-  int FalseBondDegree = 0;
-  atom *Walker = NULL;
-  atom *Candidate = NULL;
-  atom *OtherWalker = NULL;
 
   if (BondCount != 0) {
-    NoCyclicBonds = 0;
-    *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
+    *out << Verbose(1) << "Correcting Bond degree of each bond ... " << endl;
     do {
-      No = 0; // No acts as breakup flag (if 1 we still continue)
-      Walker = start;
-      while (Walker->next != end) { // go through every atom
-        Walker = Walker->next;
-
-        NoBonds = CountAtomsBonds(Walker->nr);
-        *out << Verbose(3) << "Walker " << *Walker << ": " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-        if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
-          Candidate = NULL;
-          CandidateBondNo = -1;
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
-            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-            NoBonds = CountAtomsBonds(OtherWalker->nr);
-            *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-            if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
-              if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
-                Candidate = OtherWalker;
-                CandidateBondNo = i;
-                *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
-              }
-            }
-          }
-          if ((Candidate != NULL) && (CandidateBondNo != -1)) {
-            ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
-            *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
-          } else {
-            *out << Verbose(2) << "Could not find correct degree for atom " << *Walker << "." << endl;
-            FalseBondDegree++;
-          }
-        }
-      }
+      No = SumPerAtom( &atom::CorrectBondDegree, out );
     } while (No);
     *out << " done." << endl;
@@ -245,7 +190,7 @@
     *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
   }
-  *out << FalseBondDegree << " bonds could not be corrected." << endl;
-
-  return (FalseBondDegree);
+  *out << No << " bonds could not be corrected." << endl;
+
+  return (No);
 };
 
@@ -257,5 +202,5 @@
 int molecule::CountCyclicBonds(ofstream *out)
 {
-  int No = 0;
+  NoCyclicBonds = 0;
   int *MinimumRingSize = NULL;
   MoleculeLeafClass *Subgraphs = NULL;
@@ -275,8 +220,8 @@
     Binder = Binder->next;
     if (Binder->Cyclic)
-      No++;
+      NoCyclicBonds++;
   }
   delete(BackEdgeStack);
-  return No;
+  return NoCyclicBonds;
 };
 
@@ -459,5 +404,22 @@
   }
   // set cyclic bond criterium on "same LP" basis
-  Binder = first;
+  CyclicBondAnalysis();
+
+  OutputGraphInfoPerAtom(out);
+
+  OutputGraphInfoPerBond(out);
+
+  // free all and exit
+  delete(AtomStack);
+  *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
+  return SubGraphs;
+};
+
+/** Scans through all bonds and set bond::Cyclic to true where atom::LowpointNr of both ends is equal: LP criterion.
+ */
+void molecule::CyclicBondAnalysis()
+{
+  NoCyclicBonds = 0;
+  bond *Binder = first;
   while(Binder->next != last) {
     Binder = Binder->next;
@@ -467,34 +429,34 @@
     }
   }
-
-
+};
+
+/** Output graph information per atom.
+ * \param *out output stream
+ */
+void molecule::OutputGraphInfoPerAtom(ofstream *out)
+{
   *out << Verbose(1) << "Final graph info for each atom is:" << endl;
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    *out << Verbose(2) << "Atom " << Walker->Name << " is " << ((Walker->SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
-    OutputComponentNumber(out, Walker);
-    *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
-  }
-
+  ActOnAllAtoms( &atom::OutputGraphInfo, out );
+};
+
+/** Output graph information per bond.
+ * \param *out output stream
+ */
+void molecule::OutputGraphInfoPerBond(ofstream *out)
+{
   *out << Verbose(1) << "Final graph info for each bond is:" << endl;
-  Binder = first;
+  bond *Binder = first;
   while(Binder->next != last) {
     Binder = Binder->next;
     *out << Verbose(2) << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
     *out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
-    OutputComponentNumber(out, Binder->leftatom);
+    Binder->leftatom->OutputComponentNumber(out);
     *out << " ===  ";
     *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
-    OutputComponentNumber(out, Binder->rightatom);
+    Binder->rightatom->OutputComponentNumber(out);
     *out << ">." << endl;
     if (Binder->Cyclic) // cyclic ??
       *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
   }
-
-  // free all and exit
-  delete(AtomStack);
-  *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
-  return SubGraphs;
 };
 
@@ -547,8 +509,7 @@
       Walker = BFSStack->PopFirst();
       *out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-        Binder = ListOfBondsPerAtom[Walker->nr][i];
-        if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
-          OtherAtom = Binder->GetOtherAtom(Walker);
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->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);
 #ifdef ADDHYDROGEN
           if (OtherAtom->type->Z != 1) {
@@ -662,8 +623,7 @@
           Walker = BFSStack->PopFirst();
           //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-            Binder = ListOfBondsPerAtom[Walker->nr][i];
-            if ((Binder != BackEdge) || (NumberOfBondsPerAtom[Walker->nr] == 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 = Binder->GetOtherAtom(Walker);
+          for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+            if (((*Runner) != BackEdge) || (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
+              OtherAtom = (*Runner)->GetOtherAtom(Walker);
               //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
               if (ColorList[OtherAtom->nr] == white) {
@@ -717,7 +677,7 @@
 void molecule::SetNextComponentNumber(atom *vertex, int nr)
 {
-  int i=0;
+  size_t i=0;
   if (vertex != NULL) {
-    for(;i<NumberOfBondsPerAtom[vertex->nr];i++) {
+    for(;i<vertex->ListOfBonds.size();i++) {
       if (vertex->ComponentNr[i] == -1) {   // check if not yet used
         vertex->ComponentNr[i] = nr;
@@ -727,17 +687,8 @@
         break;  // breaking here will not cause error!
     }
-    if (i == NumberOfBondsPerAtom[vertex->nr])
+    if (i == vertex->ListOfBonds.size())
       cerr << "Error: All Component entries are already occupied!" << endl;
   } else
       cerr << "Error: Given vertex is NULL!" << endl;
-};
-
-/** Output a list of flags, stating whether the bond was visited or not.
- * \param *out output stream for debugging
- */
-void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
-{
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
-    *out << vertex->ComponentNr[i] << "  ";
 };
 
@@ -751,6 +702,6 @@
     if (Walker->ComponentNr != NULL)
       Free(&Walker->ComponentNr);
-    Walker->ComponentNr = Malloc<int>(NumberOfBondsPerAtom[Walker->nr], "molecule::InitComponentNumbers: *Walker->ComponentNr");
-    for (int i=NumberOfBondsPerAtom[Walker->nr];i--;)
+    Walker->ComponentNr = Malloc<int>(Walker->ListOfBonds.size()+1, "molecule::InitComponentNumbers: *Walker->ComponentNr");
+    for (int i=Walker->ListOfBonds.size()+1;i--;)
       Walker->ComponentNr[i] = -1;
   }
@@ -763,7 +714,7 @@
 bond * molecule::FindNextUnused(atom *vertex)
 {
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
-    if (ListOfBondsPerAtom[vertex->nr][i]->IsUsed() == white)
-      return(ListOfBondsPerAtom[vertex->nr][i]);
+  for (BondList::const_iterator Runner = vertex->ListOfBonds.begin(); Runner != vertex->ListOfBonds.end(); (++Runner))
+    if ((*Runner)->IsUsed() == white)
+      return((*Runner));
   return NULL;
 };
@@ -813,5 +764,4 @@
 {
   ofstream AdjacencyFile;
-  atom *Walker = NULL;
   stringstream line;
   bool status = true;
@@ -821,12 +771,5 @@
   *out << Verbose(1) << "Saving adjacency list ... ";
   if (AdjacencyFile != NULL) {
-    Walker = start;
-    while(Walker->next != end) {
-      Walker = Walker->next;
-      AdjacencyFile << Walker->nr << "\t";
-      for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
-        AdjacencyFile << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
-      AdjacencyFile << endl;
-    }
+    ActOnAllAtoms( &atom::OutputAdjacency, &AdjacencyFile );
     AdjacencyFile.close();
     *out << Verbose(1) << "done." << endl;
@@ -850,4 +793,5 @@
   stringstream filename;
   bool status = true;
+  atom *Walker = NULL;
   char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
 
@@ -859,5 +803,5 @@
     int NonMatchNumber = 0;   // will number of atoms with differing bond structure
     int *CurrentBonds = Malloc<int>(8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
-    int CurrentBondsOfAtom;
+    size_t CurrentBondsOfAtom;
 
     // Parse the file line by line and count the bonds
@@ -871,12 +815,13 @@
       // parse into structure
       if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+        Walker = ListOfAtoms[AtomNr];
         while (!line.eof())
           line >> CurrentBonds[ ++CurrentBondsOfAtom ];
         // compare against present bonds
         //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
-        if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
-          for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
-            int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
-            int j = 0;
+        if (CurrentBondsOfAtom == Walker->ListOfBonds.size()) {
+          for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+            int id = (*Runner)->GetOtherAtom(Walker)->nr;
+            size_t j = 0;
             for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
             if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
@@ -891,5 +836,5 @@
           //out << endl;
         } else {
-          *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
+          *out << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << Walker->ListOfBonds.size() << "." << endl;
           status = false;
         }
@@ -937,9 +882,9 @@
     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
     if (Walker != NULL) // if this Walker exists in the subgraph ...
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
-        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-        if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
-          LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
-          *out << Verbose(3) << "Found local edge " << *(ListOfBondsPerAtom[Walker->nr][i]) << "." << endl;
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->nr]) { // found the bond
+          LocalStack->Push((*Runner));
+          *out << Verbose(3) << "Found local edge " << *(*Runner) << "." << endl;
           break;
         }
@@ -973,5 +918,4 @@
   class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
   atom *Walker = NULL, *OtherAtom = NULL;
-  bond *Binder = NULL;
 
   // add Root if not done yet
@@ -999,32 +943,27 @@
     // followed by n+1 till top of stack.
     Walker = AtomStack->PopFirst(); // pop oldest added
-    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
-    for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-      Binder = ListOfBondsPerAtom[Walker->nr][i];
-      if (Binder != NULL) { // don't look at bond equal NULL
-        OtherAtom = Binder->GetOtherAtom(Walker);
-        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
+    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << Walker->ListOfBonds.size() << " bonds." << endl;
+    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      if ((*Runner) != NULL) { // don't look at bond equal NULL
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *(*Runner) << "." << endl;
         if (ColorList[OtherAtom->nr] == white) {
-          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
+          if ((*Runner) != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
             ColorList[OtherAtom->nr] = lightgray;
           PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
           ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
           *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
-          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
+          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && ((*Runner) != Bond))) ) { // Check for maximum distance
             *out << Verbose(3);
             if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
               AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
               *out << "Added OtherAtom " << OtherAtom->Name;
-              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-              AddedBondList[Binder->nr]->Type = Binder->Type;
-              *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
+              AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
+              *out << " and bond " << *(AddedBondList[(*Runner)->nr]) << ", ";
             } else {  // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
               *out << "Not adding OtherAtom " << OtherAtom->Name;
-              if (AddedBondList[Binder->nr] == NULL) {
-                AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-                AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-                AddedBondList[Binder->nr]->Type = Binder->Type;
-                *out << ", added Bond " << *(AddedBondList[Binder->nr]);
+              if (AddedBondList[(*Runner)->nr] == NULL) {
+                AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
+                *out << ", added Bond " << *(AddedBondList[(*Runner)->nr]);
               } else
                 *out << ", not added Bond ";
@@ -1033,20 +972,18 @@
             AtomStack->Push(OtherAtom);
           } else { // out of bond order, then replace
-            if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
+            if ((AddedAtomList[OtherAtom->nr] == NULL) && ((*Runner)->Cyclic))
               ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
-            if (Binder == Bond)
+            if ((*Runner) == Bond)
               *out << Verbose(3) << "Not Queueing, is the Root bond";
             else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
               *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
-            if (!Binder->Cyclic)
+            if (!(*Runner)->Cyclic)
               *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
-            if (AddedBondList[Binder->nr] == NULL) {
+            if (AddedBondList[(*Runner)->nr] == NULL) {
               if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
-                AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-                AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-                AddedBondList[Binder->nr]->Type = Binder->Type;
+                AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
               } else {
 #ifdef ADDHYDROGEN
-                if (!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
+                if (!Mol->AddHydrogenReplacementAtom(out, (*Runner), AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
                   exit(1);
 #endif
@@ -1057,12 +994,10 @@
           *out << Verbose(3) << "Not Adding, has already been visited." << endl;
           // This has to be a cyclic bond, check whether it's present ...
-          if (AddedBondList[Binder->nr] == NULL) {
-            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
-              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-              AddedBondList[Binder->nr]->Type = Binder->Type;
+          if (AddedBondList[(*Runner)->nr] == NULL) {
+            if (((*Runner) != Bond) && ((*Runner)->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
+              AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
             } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
 #ifdef ADDHYDROGEN
-              if(!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
+              if(!Mol->AddHydrogenReplacementAtom(out, (*Runner), AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
                 exit(1);
 #endif
@@ -1080,4 +1015,19 @@
   delete(AtomStack);
 };
+
+/** Adds a bond as a copy to a given one
+ * \param *left leftatom of new bond
+ * \param *right rightatom of new bond
+ * \param *CopyBond rest of fields in bond are copied from this
+ * \return pointer to new bond
+ */
+bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
+{
+  bond *Binder = AddBond(left, right, CopyBond->BondDegree);
+  Binder->Cyclic = CopyBond->Cyclic;
+  Binder->Type = CopyBond->Type;
+  return Binder;
+};
+
 
 /** Adds bond structure to this molecule from \a Father molecule.
@@ -1122,9 +1072,9 @@
         status = false;
       } else {
-        for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
-          OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          OtherAtom = (*Runner)->GetOtherAtom(Walker);
           if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
-            *out << Verbose(4) << "Endpoints of Bond " << Father->ListOfBondsPerAtom[Walker->nr][i] << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
-            AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], Father->ListOfBondsPerAtom[Walker->nr][i]->BondDegree);
+            *out << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
+            AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
           }
         }
@@ -1164,6 +1114,6 @@
       for(KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
         Walker2 = FindAtom(*runners);
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          if (ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker) == Walker2) {
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner)->GetOtherAtom(Walker) == Walker2) {
             BondStatus = true;
             break;
Index: src/molecule_template.hpp
===================================================================
--- src/molecule_template.hpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/molecule_template.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -153,4 +153,90 @@
 };
 
+// ========================= Summing over each Atoms =================================== //
+
+// zero arguments
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() const )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() const ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+// one argument
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T), T t )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T) const, T t )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T), T t ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T) const, T t ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+
+
 // ================== Acting with each Atoms on same molecule ========================== //
 
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/moleculelist.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -376,4 +376,5 @@
   atom *Walker = NULL;
   atom *Runner = NULL;
+  bond *Binder = NULL;
   double ***FitConstant = NULL, **correction = NULL;
   int a, b;
@@ -482,5 +483,5 @@
     while (Walker->next != (*ListRunner)->end) {
       Walker = Walker->next;
-      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      //cout << 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
@@ -488,7 +489,8 @@
         while (Runner->next != (*ListRunner)->end) {
           Runner = Runner->next;
-          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
-          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one 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!)
             // 4. evaluate the morse potential for each matrix component and add up
             distance = Runner->x.Distance(&Walker->x);
@@ -812,6 +814,7 @@
 bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
 {
-  atom *Walker = NULL, *OtherWalker = NULL;
-  bond *Binder = NULL;
+  atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  atom *Father = NULL;
   bool status = true;
   int AtomNo;
@@ -830,18 +833,17 @@
     while (Walker->next != Leaf->end) {
       Walker = Walker->next;
-      AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker
-      for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all
-        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
-        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
+      Father = Walker->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
         if (OtherWalker != NULL) {
           if (OtherWalker->nr > Walker->nr)
-            Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
+            Leaf->AddBond(Walker, OtherWalker, (*Runner)->BondDegree);
         } else {
-          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
+          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
           status = false;
         }
       }
     }
-    Leaf->CreateListOfBondsPerAtom(out);
     FragmentCounter++;
     if (next != NULL)
Index: src/unittests/Makefile.am
===================================================================
--- src/unittests/Makefile.am	(revision b8b75d6a56dc001347024c088c6590bec244d5c7)
+++ src/unittests/Makefile.am	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -1,7 +1,7 @@
 INCLUDES = -I$(top_srcdir)/src
 
-noinst_PROGRAMS =  ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest VectorUnitTest
+noinst_PROGRAMS =  ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest ListOfBondsUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest VectorUnitTest
 
-TESTS = ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest VectorUnitTest
+TESTS = ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest ListOfBondsUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest VectorUnitTest
 check_PROGRAMS = $(TESTS)
 
@@ -26,4 +26,9 @@
 AnalysisPairCorrelationUnitTest_LDADD = ../libmolecuilder.a
 
+ListOfBondsUnitTest_SOURCES = listofbondsunittest.cpp listofbondsunittest.hpp
+ListOfBondsUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
+ListOfBondsUnitTest_LDFLAGS = $(CPPUNIT_LIBS) -ldl
+ListOfBondsUnitTest_LDADD = ../libmolecuilder.a
+
 VectorUnitTest_SOURCES = defs.hpp helpers.hpp leastsquaremin.hpp memoryallocator.hpp memoryusageobserver.hpp vectorunittest.cpp vectorunittest.hpp vector.hpp verbose.hpp 
 VectorUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
Index: src/unittests/listofbondsunittest.cpp
===================================================================
--- src/unittests/listofbondsunittest.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
+++ src/unittests/listofbondsunittest.cpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -0,0 +1,270 @@
+/*
+ * listofbondsunittest.cpp
+ *
+ *  Created on: 18 Oct 2009
+ *      Author: user
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "listofbondsunittest.hpp"
+
+#include "atom.hpp"
+#include "bond.hpp"
+#include "element.hpp"
+#include "molecule.hpp"
+#include "periodentafel.hpp"
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( ListOfBondsTest );
+
+
+void ListOfBondsTest::setUp()
+{
+  atom *Walker = NULL;
+
+  // init private all pointers to zero
+  TestMolecule = NULL;
+  hydrogen = NULL;
+  tafel = NULL;
+
+  // construct element
+  hydrogen = new element;
+  hydrogen->Z = 1;
+  strcpy(hydrogen->name, "hydrogen");
+  strcpy(hydrogen->symbol, "H");
+
+
+  // construct periodentafel
+  tafel = new periodentafel;
+  tafel->AddElement(hydrogen);
+
+  // construct molecule (tetraeder of hydrogens)
+  TestMolecule = new molecule(tafel);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(1., 0., 1. );
+  TestMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(0., 1., 1. );
+  TestMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(1., 1., 0. );
+  TestMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(0., 0., 0. );
+  TestMolecule->AddAtom(Walker);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
+
+};
+
+
+void ListOfBondsTest::tearDown()
+{
+  // remove
+  delete(TestMolecule);
+  // note that all the atoms are cleaned by TestMolecule
+  delete(tafel);
+  // note that element is cleaned by periodentafel
+};
+
+/** Unit Test of molecule::AddBond()
+ *
+ */
+void ListOfBondsTest::AddingBondTest()
+{
+  bond *Binder = NULL;
+  atom *atom1 = TestMolecule->start->next;
+  atom *atom2 = atom1->next;
+  CPPUNIT_ASSERT( atom1 != NULL );
+  CPPUNIT_ASSERT( atom2 != NULL );
+
+  // add bond
+  Binder = TestMolecule->AddBond(atom1, atom2, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+  bond *TestBond = TestMolecule->first->next;
+  CPPUNIT_ASSERT_EQUAL ( TestBond, Binder );
+
+  // check that bond contains the two atoms
+  CPPUNIT_ASSERT_EQUAL( true, Binder->Contains(atom1) );
+  CPPUNIT_ASSERT_EQUAL( true, Binder->Contains(atom2) );
+
+  // check that bond is present in both atoms
+  bond *TestBond1 = *(atom1->ListOfBonds.begin());
+  CPPUNIT_ASSERT_EQUAL( TestBond1, Binder );
+  bond *TestBond2 = *(atom2->ListOfBonds.begin());
+  CPPUNIT_ASSERT_EQUAL( TestBond2, Binder );
+};
+
+/** Unit Test of molecule::RemoveBond()
+ *
+ */
+void ListOfBondsTest::RemovingBondTest()
+{
+  bond *Binder = NULL;
+  atom *atom1 = TestMolecule->start->next;
+  atom *atom2 = atom1->next;
+  CPPUNIT_ASSERT( atom1 != NULL );
+  CPPUNIT_ASSERT( atom2 != NULL );
+
+  // add bond
+  Binder = TestMolecule->AddBond(atom1, atom2, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+
+  // remove bond
+  TestMolecule->RemoveBond(Binder);
+
+  // check if removed from atoms
+  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom2->ListOfBonds.size() );
+
+  // check if removed from molecule
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->first->next, TestMolecule->last );
+};
+
+/** Unit Test of molecule::RemoveBonds()
+ *
+ */
+void ListOfBondsTest::RemovingBondsTest()
+{
+  bond *Binder = NULL;
+  atom *atom1 = TestMolecule->start->next;
+  atom *atom2 = atom1->next;
+  atom *atom3 = atom2->next;
+  CPPUNIT_ASSERT( atom1 != NULL );
+  CPPUNIT_ASSERT( atom2 != NULL );
+  CPPUNIT_ASSERT( atom3 != NULL );
+
+  // add bond
+  Binder = TestMolecule->AddBond(atom1, atom2, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+  Binder = TestMolecule->AddBond(atom1, atom3, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+  Binder = TestMolecule->AddBond(atom2, atom3, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+
+  // 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() );
+
+  // remove bond
+  TestMolecule->RemoveBonds(atom1);
+
+  // 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() );
+
+  // check if removed from molecule
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->first->next, Binder );
+  CPPUNIT_ASSERT_EQUAL( Binder->next, TestMolecule->last );
+};
+
+/** Unit Test of delete(bond *)
+ *
+ */
+void ListOfBondsTest::DeleteBondTest()
+{
+  bond *Binder = NULL;
+  atom *atom1 = TestMolecule->start->next;
+  atom *atom2 = atom1->next;
+  CPPUNIT_ASSERT( atom1 != NULL );
+  CPPUNIT_ASSERT( atom2 != NULL );
+
+  // add bond
+  Binder = TestMolecule->AddBond(atom1, atom2, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+
+  // remove bond
+  delete(Binder);
+
+  // check if removed from atoms
+  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom2->ListOfBonds.size() );
+
+  // check if removed from molecule
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->first->next, TestMolecule->last );
+};
+
+/** Unit Test of molecule::RemoveAtom()
+ *
+ */
+void ListOfBondsTest::RemoveAtomTest()
+{
+  bond *Binder = NULL;
+  atom *atom1 = TestMolecule->start->next;
+  atom *atom2 = atom1->next;
+  CPPUNIT_ASSERT( atom1 != NULL );
+  CPPUNIT_ASSERT( atom2 != NULL );
+
+  // add bond
+  Binder = TestMolecule->AddBond(atom1, atom2, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+
+  // remove atom2
+  TestMolecule->RemoveAtom(atom2);
+
+  // check bond if removed from other atom
+  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
+
+  // check if removed from molecule
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->first->next, TestMolecule->last );
+};
+
+/** Unit Test of delete(atom *)
+ *
+ */
+void ListOfBondsTest::DeleteAtomTest()
+{
+  bond *Binder = NULL;
+  atom *atom1 = TestMolecule->start->next;
+  atom *atom2 = atom1->next;
+  CPPUNIT_ASSERT( atom1 != NULL );
+  CPPUNIT_ASSERT( atom2 != NULL );
+
+  // add bond
+  Binder = TestMolecule->AddBond(atom1, atom2, 1);
+  CPPUNIT_ASSERT( Binder != NULL );
+
+  // remove atom2
+  delete(atom2);
+
+  // check bond if removed from other atom
+  CPPUNIT_ASSERT_EQUAL( (size_t) 0, atom1->ListOfBonds.size() );
+
+  // check if removed from molecule
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->first->next, TestMolecule->last );
+};
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/listofbondsunittest.hpp
===================================================================
--- src/unittests/listofbondsunittest.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
+++ src/unittests/listofbondsunittest.hpp	(revision 266237a46efdbba0be37a76f5e86b6814aac18d7)
@@ -0,0 +1,48 @@
+/*
+ * listofbondsunittest.hpp
+ *
+ *  Created on: 18 Oct 2009
+ *      Author: user
+ */
+
+#ifndef LISTOFBONDSUNITTEST_HPP_
+#define LISTOFBONDSUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class element;
+class molecule;
+class periodentafel;
+
+/********************************************** Test classes **************************************/
+
+class ListOfBondsTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( ListOfBondsTest) ;
+    CPPUNIT_TEST ( AddingBondTest );
+    CPPUNIT_TEST ( RemovingBondTest );
+    CPPUNIT_TEST ( RemovingBondsTest );
+    CPPUNIT_TEST ( RemoveAtomTest );
+    CPPUNIT_TEST ( DeleteBondTest );
+    CPPUNIT_TEST ( DeleteAtomTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void AddingBondTest();
+      void RemovingBondTest();
+      void RemovingBondsTest();
+      void RemoveAtomTest();
+      void DeleteBondTest();
+      void DeleteAtomTest();
+
+private:
+
+      molecule *TestMolecule;
+      element *hydrogen;
+      periodentafel *tafel;
+};
+
+
+#endif /* LISTOFBONDSUNITTEST_HPP_ */
