Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 683914b4876eed83b28d109437ddb9f72c5a18f5)
+++ src/molecules.cpp	(revision 41baafea7091ff08639dc60305b81c5e1871ee8e)
@@ -1442,6 +1442,6 @@
 void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
 {
-  atom *Walker = NULL, *OtherWalker = NULL;
-  int No, NoBonds;
+  atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
+  int No, NoBonds, CandidateBondNo;
   int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
   molecule **CellList;
@@ -1556,7 +1556,8 @@
     CreateListOfBondsPerAtom(out);
   		
-    // correct Bond degree of each bond by checking of updated(!) sum of bond degree for an atom match its valence count
-    // bond degrres are correctled iteratively by one, so that 2-2 instead of 1-3 or 3-1 corrections are favoured: We want
-    // a rather symmetric distribution of higher bond degrees
+    // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
+    // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
+    // preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
+    // double bonds as was expected.
 		if (BondCount != 0) {
       NoCyclicBonds = 0;
@@ -1567,20 +1568,30 @@
         while (Walker->next != end) { // go through every atom
           Walker = Walker->next;
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
-            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-	          // count valence of first partner (updated!), might have changed during last bond partner
-	          NoBonds = 0;
-            for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
-              NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-	          *out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-	          if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check NoBonds of other atom
+          // count valence of first partner
+          NoBonds = 0;
+          for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
+            NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
+          *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;
+            for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
+              OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
 	            // count valence of second partner
               NoBonds = 0;
               for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
                 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
-              *out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-	            if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) // increase bond degree by one
-                ListOfBondsPerAtom[Walker->nr][i]->BondDegree++;
+              *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) {
+              ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
+              *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
+            }
           }
 		    }
@@ -1801,4 +1812,8 @@
 
 /** Analyses the cycles found and returns minimum of all cycle lengths.
+ * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root,
+ * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
+ * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
+ * as cyclic and print out the cycles. 
  * \param *out output stream for debugging
  * \param *BackEdgeStack stack with all back edges found during DFS scan
@@ -1812,5 +1827,5 @@
   enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
   class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring 
-  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms 
+  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
   bond *Binder = NULL, *BackEdge = NULL;
@@ -1843,36 +1858,72 @@
     BFSStack->Push(Walker);
     TouchedStack->Push(Walker);
-    //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
+    *out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
     OtherAtom = NULL;
-    while ((Walker != Root) && ((OtherAtom == NULL) || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]))) {  // look for Root
+    do {  // look for Root
       Walker = BFSStack->PopFirst();
-      //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
+      *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);
-          //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
-          if (ColorList[OtherAtom->nr] == white) {
-            TouchedStack->Push(OtherAtom);
-            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 << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
-            if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
-              //*out << Verbose(3) << "Putting OtherAtom into queue." << endl;
-              BFSStack->Push(OtherAtom);
+#ifdef ADDHYDROGEN          
+          if (OtherAtom->type->Z != 1) {
+#endif
+            *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
+            if (ColorList[OtherAtom->nr] == white) {
+              TouchedStack->Push(OtherAtom);
+              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 << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+              //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
+                *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
+                BFSStack->Push(OtherAtom);
+              //}
+            } else {
+              *out << Verbose(3) << "Not Adding, has already been visited." << endl;
             }
+            if (OtherAtom == Root)
+              break;
+#ifdef ADDHYDROGEN          
           } else {
-            //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
+            *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
+            ColorList[OtherAtom->nr] = black;
           }
+#endif
         } else {
-          //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
+          *out << Verbose(2) << "Bond " << *Binder << " not Visiting, is the back edge." << endl;
         }
       }
       ColorList[Walker->nr] = black;
-      //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
-    }
+      *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
+      if (OtherAtom == Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
+        // step through predecessor list
+        while (OtherAtom != BackEdge->rightatom) {
+          if (!OtherAtom->GetTrueFather()->IsCyclic)  // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
+            break;
+          else
+            OtherAtom = PredecessorList[OtherAtom->nr];
+        }
+        if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
+          *out << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl;\
+          do {
+            OtherAtom = TouchedStack->PopLast();
+            if (PredecessorList[OtherAtom->nr] == Walker) {
+              *out << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl;
+              PredecessorList[OtherAtom->nr] = NULL;
+              ShortestPathList[OtherAtom->nr] = -1;
+              ColorList[OtherAtom->nr] = white;
+              BFSStack->RemoveItem(OtherAtom);
+            }
+          } while ((!TouchedStack->IsEmpty()) && (PredecessorList[OtherAtom->nr] == NULL));
+          TouchedStack->Push(OtherAtom);  // last was wrongly popped
+          OtherAtom = BackEdge->rightatom; // set to not Root
+        } else
+          OtherAtom = Root;
+      }
+    } while ((OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
     
-    if (Walker == Root) {
+    if (OtherAtom == Root) {
       // now climb back the predecessor list and thus find the cycle members
       NumCycles++;
@@ -1880,4 +1931,5 @@
       Root->GetTrueFather()->IsCyclic = true;
       *out << Verbose(1) << "Found ring contains: ";
+      Walker = Root;
       while (Walker != BackEdge->rightatom) {
         *out << Walker->Name << " <-> ";
