Index: c/border.cpp
===================================================================
--- src/border.cpp	(revision f683feb0f2b4232b7bc42f3893078413c798f659)
+++ 	(revision )
@@ -1,302 +1,0 @@
-#include "molecules.hpp"
-#include "boundary.hpp"
-
-
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol)
-{
-  /* d2 ist der Normalenvektor auf dem Dreieck,
-   * d1 ist der Vektor, der normal auf der Kante und d2 steht.
-   */
-  Vector dif_a;
-  Vector dif_b;
-  Vector Chord;
-  Vector AngleCheck;
-  atom *Walker;
-
-  dif_a.CopyVector(&a.x);
-  dif_a.SubtractVector(&Candidate.x);
-  dif_b.CopyVector(&b.x);
-  dif_b.SubtractVector(&Candidate.x);
-  Chord.CopyVector(&a.x);
-  Chord.SubtractVector(&b.x);
-  AngleCheck.CopyVector(&dif_a);
-  AngleCheck.ProjectOntoPlane(&Chord);
-
-  //Storage eintrag fuer aktuelles Atom
-  if (Chord.Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
-  {
-
-    if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
-      {
-	Storage[0]=(double)Candidate.nr;
-	Storage[1]=1;
-	Storage[2]=AngleCheck.Angle(d2);
-      }
-    else
-      {
-	if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
-	    (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
-	  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-	  {
-	    Storage[0]=(double)Candidate.nr;
-	    Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-	    Storage[2]=AngleCheck.Angle(d2);
-	  }
-      }
-  }
-
-  if (n<5)
-    {
-      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
-	{
-	  while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
-	    {
-	      Walker = Walker->next;
-	    }
-
-	  Find_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
-	}
-    }
-};
-
-
-void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
-{
-  Vector CenterOfLine = Line.endpoints[0]->node->x;
-  Vector direction1;
-  Vector direction2;
-  Vector helper;
-  atom* Walker;
-
-  double Storage[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
-
-  
-  helper.CopyVector(&(Line.endpoints[0]->node->x));
-  for (int i =0; i<3; i++)
-    {
-      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr!=Line.endpoints[1]->node->nr)
-	{
-	  helper.SubtractVector(&T.endpoints[i]->node->x);
-	  break;
-	}
-    }
-
-
-  direction1.CopyVector(&Line.endpoints[0]->node->x);
-  direction1.SubtractVector(&Line.endpoints[1]->node->x);
-  direction1.VectorProduct(T.NormalVector);
-
-  if (direction1.ScalarProduct(&helper)>0)
-    {
-      direction1.Scale(-1);
-    }
-
-  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &direction1, T.NormalVector, Storage, RADIUS, *mol);
-  
-  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
-  // Next Triangle is Line, atom with number in Storage[0]
-
-  Walker= mol->start;
-  while (Walker->nr != (int)Storage[0])
-    {
-      Walker = Walker->next;
-    }
-
-  AddPoint(Walker);
-
-  BPS[0] = new class BoundaryPointSet(Walker);
-  BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
-  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(Walker);
-  BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
-  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BLS[2] = new class BoundaryLineSet(Line);
-
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
-  TrianglesOnBoundaryCount++;
-
-  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
-    {
-      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
-	{
-	  LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
-	  LinesOnBoundaryCount++;
-	}
-    }
-  BTS->GetNormalVector(*BTS->NormalVector);
-
-  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
-      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
-    {
-      BTS->NormalVector->Scale(-1);
-    }
-  
-};
-
-
-void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
-{
-  int i;
-  Vector *AngleCheck;
-  atom* Walker;
-
-  AngleCheck->CopyVector(&Candidate.x);
-  AngleCheck->SubtractVector(&a.x);
-  if (AngleCheck->ScalarProduct(&Oben) < Storage[1])
-    {
-      Storage[0]=(double)(Candidate.nr);
-      Storage[1]=AngleCheck->ScalarProduct(&Oben);
-    };
-
-  if (n<5)
-    {
-      for (i = 0; i< mol.NumberOfBondsPerAtom[Candidate.nr]; i++)
-	{
-	  Walker = mol.start;
-	  while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
-	    {
-	      Walker = Walker->next;
-	    };
-	  
-	  Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
-	    };
-    };
-  
-
-};
-
-
-void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
-{
-  int i=0;
-  atom Walker;
-  atom Walker2;
-  atom Walker3;
-  int max_index[3];
-  double max_coordinate[3];
-  Vector Oben;
-  Vector helper;
-
-  Oben.Zero();
-
-
-  for(i =0; i<3; i++)
-    {
-      max_index[i] =-1;
-      max_coordinate[i] =-1;
-    }
-
-  Walker = *mol.start;
-  while (Walker.next != NULL)
-    {
-      for (i=0; i<3; i++)
-	{
-	  if (Walker.x.x[i]>max_coordinate[i])
-	    {
-	      max_coordinate[i]=Walker.x.x[i];
-	      max_index[i]=Walker.nr;
-	    }
-	}
-    }
-
-  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-  const int k=0;
-
-  Oben.x[k]=1;
-  Walker = *mol.start;
-  while (Walker.nr != max_index[k])
-    {
-      Walker = *Walker.next;
-    }
-
-  double Storage[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[2]=-10.;  // This will be an angle looking for the third point.
-
-
-  for (i=0; i< mol.NumberOfBondsPerAtom[Walker.nr]; i++)
-    {
-      Walker2 = *mol.start;
-      while (Walker2.nr != (mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom->nr : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr) ) 
-	{
-	  Walker2 = *Walker2.next;
-	}
-
-      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
-    }
-
-  Walker2 = *mol.start;
-
-  while (Walker2.nr != int(Storage[0]))
-    {
-      Walker = *Walker.next;
-    }
-  
-  helper.CopyVector(&Walker.x);
-  helper.SubtractVector(&Walker2.x);
-  Oben.ProjectOntoPlane(&helper);
-  helper.VectorProduct(&Oben);
-
-       Find_next_suitable_point(Walker, Walker2, *(mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom), 0, &helper, &Oben, Storage, RADIUS, mol);
-  Walker3 = *mol.start;
-  while (Walker3.nr != int(Storage[0]))
-    {
-      Walker3 = *Walker3.next;
-    }
-
-  //Starting Triangle is Walker, Walker2, index Storage[0]
-
-  AddPoint(&Walker);
-  AddPoint(&Walker2);
-  AddPoint(&Walker3);
-
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker2);
-  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker3);
-  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker2);
-  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
-
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
-  TrianglesOnBoundaryCount++;
-
-  for(int i=0;i<NDIM;i++) 
-    {
-      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
-      LinesOnBoundaryCount++;
-    };
-
-       BTS->GetNormalVector(*BTS->NormalVector);
-
-       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
-	 {
-	   BTS->NormalVector->Scale(-1);
-	 }
-};
-
-
-void Find_non_convex_border(Tesselation* Tess, molecule mol) 
-{
-  const double RADIUS =6;
-  Tess->Find_starting_triangle(mol, RADIUS);
-
-  for (LineMap::iterator baseline = Tess->LinesOnBoundary.begin(); baseline != Tess->LinesOnBoundary.end(); baseline++) 
-    if (baseline->second->TrianglesCount == 1) 
-      {
-	Tess->Find_next_suitable_triangle(&mol, baseline->second, baseline->second->triangles.begin()->second, RADIUS); //the line is there, so there is a triangle, but only one.
-
-      }
-    else
-      {
-	printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
-      }
-};
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision f683feb0f2b4232b7bc42f3893078413c798f659)
+++ src/boundary.cpp	(revision a8bcea6105d46c2f1e551ab69003f9d1cc03e44f)
@@ -456,4 +456,40 @@
 
 
+/*
+ * This function creates the tecplot file, displaying the tesselation of the hull.
+ * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ */
+void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol)
+{
+  // 8. Store triangles in tecplot file
+  if (tecplot != NULL) {
+    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+    *tecplot << "ZONE T=\"TRIANGLES\", N=" <<  TesselStruct->PointsOnBoundaryCount << ", E=" <<  TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+    int *LookupList = new int[mol->AtomCount];
+    for (int i=0;i<mol->AtomCount;i++)
+      LookupList[i] = -1;
+
+    // print atom coordinates
+    *out << Verbose(2) << "The following triangles were created:";
+    int Counter = 1;
+    atom *Walker = NULL;
+    for (PointMap::iterator target =  TesselStruct->PointsOnBoundary.begin(); target !=  TesselStruct->PointsOnBoundary.end(); target++) {
+      Walker = target->second->node;
+      LookupList[Walker->nr] = Counter++;
+      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+    }
+    *tecplot << endl;
+      // print connectivity
+    for (TriangleMap::iterator runner =  TesselStruct->TrianglesOnBoundary.begin(); runner !=  TesselStruct->TrianglesOnBoundary.end(); runner++) {
+      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
+      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+    }
+    delete[](LookupList);
+    *out << endl;
+  }
+}
+
 /** Determines the volume of a cluster.
  * Determines first the convex envelope, then tesselates it and calculates its volume.
@@ -478,5 +514,5 @@
   double a,b,c;
 
-  Find_non_convex_border(*TesselStruct, mol);
+  Find_non_convex_border(out, tecplot, *TesselStruct, mol);
   /*
   // 1. calculate center of gravity
@@ -560,32 +596,6 @@
 
 
-
-  // 8. Store triangles in tecplot file
-  if (tecplot != NULL) {
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES\", N=" <<  TesselStruct->PointsOnBoundaryCount << ", E=" <<  TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int *LookupList = new int[mol->AtomCount];
-    for (int i=0;i<mol->AtomCount;i++)
-      LookupList[i] = -1;
-
-    // print atom coordinates
-    *out << Verbose(2) << "The following triangles were created:";
-    int Counter = 1;
-    atom *Walker = NULL;
-    for (PointMap::iterator target =  TesselStruct->PointsOnBoundary.begin(); target !=  TesselStruct->PointsOnBoundary.end(); target++) {
-      Walker = target->second->node;
-      LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
-    }
-    *tecplot << endl;
-      // print connectivity
-    for (TriangleMap::iterator runner =  TesselStruct->TrianglesOnBoundary.begin(); runner !=  TesselStruct->TrianglesOnBoundary.end(); runner++) {
-      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
-    }
-    delete[](LookupList);
-    *out << endl;
-  }
+  write_tecplot_file(out, tecplot, TesselStruct, mol);
+
 
   // free reference lists
@@ -699,4 +709,5 @@
 Tesselation::~Tesselation()
 {
+	cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
   for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
     delete(runner->second);
@@ -1028,5 +1039,5 @@
  *  with all neighbours of the candidate.
  */
-void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule* mol)
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
 {
   /* d2 is normal vector on the triangle
@@ -1045,41 +1056,60 @@
   AngleCheck.ProjectOntoPlane(Chord);
 
-  if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  if (problem)
   {
-
-    if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
+	  cout << "Atom number" << Candidate->nr << endl;
+	  Candidate->x.Output((ofstream *)&cout);
+	  cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
+  }
+
+  if (a != Candidate and b != Candidate)
+  {
+
+	  if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+	  {
+		  if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[0]) //This will give absolute preference to those in "right-hand" quadrants
+		  {
+			  Opt_Candidate = Candidate;
+			  Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+			  Storage[1]=AngleCheck.Angle(d2);
+		  }
+		  else
+		  {
+			  if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]< AngleCheck.Angle(d2)) or \
+					  (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]> AngleCheck.Angle(d2)))
+				  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+			  {
+				  Opt_Candidate = Candidate;
+				  Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+				  Storage[1]=AngleCheck.Angle(d2);
+			  }
+			  else
+			  {
+				  if (problem)
+					  cout << "Loses to better candidate" << endl;
+			  }
+		  }
+	  }
+	  else
+	  {
+		  if (problem)
+			  cout << "erfuellt Dreiecksbedingung fuer sehne nicht" <<endl;
+	  }
+  }
+  else
+  {
+	  if (problem)
+	  		  cout << "identisch mit Ursprungslinie" << endl;
+  }
+
+  if (n<5) // Five is the recursion level threshold.
+  {
+	  for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
       {
-    	Storage[0]=(double)Candidate->nr;
-    	Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-    	Storage[2]=AngleCheck.Angle(d2);
+		  Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+
+		  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
       }
-    else
-      {
-	if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
-	    (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
-	  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-	  {
-	    Storage[0]=(double)Candidate->nr;
-	    Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-	    Storage[2]=AngleCheck.Angle(d2);
-	  }
-      }
-  }
-
-  if (n<5) // Five is the recursion level threshold.
-    {
-      for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
-	{
-		Walker= mol->start; // go through all atoms
-
-	  while (Walker->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr))
-	    { // until atom found which belongs to bond
-	      Walker = Walker->next;
-	    }
-
-
-	  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Storage, RADIUS, mol);  //call function again
-	}
-    }
+  }
 };
 
@@ -1087,7 +1117,7 @@
  * this function fins a triangle to a line, adjacent to an existing one.
  */
-void Tesselation::Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
-{
-	printf("Looking for next suitable triangle \n");
+void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS)
+{
+	cout << Verbose(1) << "Looking for next suitable triangle \n";
   Vector direction1;
   Vector helper;
@@ -1095,11 +1125,11 @@
   atom* Walker;
 
-  double *Storage;
-  Storage = new double[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
-
-
+  double Storage[2];
+  Storage[0]=-2.;    // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[1]=9999999.;  // This is also lower then any value produced by an eligible atom, which are all positive
+  atom* Opt_Candidate = NULL;
+
+
+  cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
   helper.CopyVector(&(Line.endpoints[0]->node->x));
   for (int i =0; i<3; i++)
@@ -1115,5 +1145,5 @@
   direction1.CopyVector(&Line.endpoints[0]->node->x);
   direction1.SubtractVector(&Line.endpoints[1]->node->x);
-  direction1.VectorProduct(T.NormalVector);
+  direction1.VectorProduct(&(T.NormalVector));
 
   if (direction1.ScalarProduct(&helper)>0)
@@ -1125,22 +1155,30 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-printf("Looking for third point candidates for triangle \n");
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, T.NormalVector, Storage, RADIUS, mol);
-
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
+
+  if (Opt_Candidate==NULL)
+  {
+	  cout << Verbose(1) << "No new Atom found, triangle construction will crash" << endl;
+	  write_tecplot_file(out, tecplot, this, mol);
+	  tecplot->flush();
+	  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol, 1);
+
+  }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
-  // Next Triangle is Line, atom with number in Storage[0]
-
-  Walker= mol->start;
-  while (Walker->nr != (int)Storage[0])
-    {
-      Walker = Walker->next;
-    }
-
-  AddPoint(Walker);
-
-  BPS[0] = new class BoundaryPointSet(Walker);
+
+  cout << Verbose(1) << "Adding exactly one Walker for reasons completely unknown to me ... " << endl;
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  AddPoint(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(Walker);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
@@ -1151,4 +1189,5 @@
   TrianglesOnBoundaryCount++;
 
+  cout << Verbose(1) << "Checking for already present lines ... " << endl;
   for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
     {
@@ -1160,44 +1199,40 @@
 	}
     }
-  BTS->GetNormalVector(*BTS->NormalVector);
-
-  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
-      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
+  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
+
+  BTS->GetNormalVector(BTS->NormalVector);
+
+  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector))<0 && Storage[0]>0) ||
+		  ((BTS->NormalVector.ScalarProduct(&(T.NormalVector))>0 && Storage[0]<0)) )
     {
-      BTS->NormalVector->Scale(-1);
-    }
-
-};
-
-
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, double Storage[3], molecule* mol)
-{
-	printf("Looking for second point of starting triangle \n");
-  int i;
-  Vector *AngleCheck;
-  atom* Walker;
-
-  AngleCheck->CopyVector(&(Candidate->x));
-  AngleCheck->SubtractVector(&(a->x));
-  if (AngleCheck->Angle(&Oben) < Storage[1])
+      BTS->NormalVector.Scale(-1);
+    };
+
+};
+
+
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, atom*& Opt_Candidate, double Storage[2], molecule* mol, double RADIUS)
+{
+	cout << Verbose(1) << "Looking for second point of starting triangle, recursive level "<< n <<endl;;
+	int i;
+	Vector AngleCheck;
+	atom* Walker;
+
+	AngleCheck.CopyVector(&(Candidate->x));
+	AngleCheck.SubtractVector(&(a->x));
+	if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
     {
-	  //printf("Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-      Storage[0]=(double)(Candidate->nr);
-      Storage[1]=AngleCheck->Angle(&Oben);
-      //printf("Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+	  //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+      Opt_Candidate=Candidate;
+      Storage[0]=AngleCheck.Angle(&Oben);
+      //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
     };
-  printf("%d \n", n);
-  if (n<5)
+	if (n<5)
     {
       for (i = 0; i< mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
-	{
-	  Walker = mol->start;
-	  while (Candidate->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr))
-	    {
-	      Walker = Walker->next;
-	    };
-
-	  Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Storage, mol);
-	    };
+      {
+    	  Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+    	  Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+      };
     };
 
@@ -1208,5 +1243,5 @@
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-	printf("Looking for starting triangle \n");
+	cout << Verbose(1) << "Looking for starting triangle \n";
   int i=0;
   atom* Walker;
@@ -1227,5 +1262,5 @@
       max_coordinate[i] =-1;
     }
-printf("Molecule mol is there and has %d Atoms \n", mol->AtomCount);
+cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
   Walker = mol->start;
   while (Walker->next != mol->end)
@@ -1242,5 +1277,5 @@
     }
 
-printf("Found starting atom \n");
+cout << Verbose(1) << "Found starting atom \n";
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
   const int k=2;
@@ -1253,54 +1288,40 @@
       Walker = Walker->next;
     }
-printf("%d \n", Walker->nr);
-  double Storage[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[2]=999999.;  // This will be an angle looking for the third point.
-printf("%d \n", mol->NumberOfBondsPerAtom[Walker->nr]);
+cout << Verbose(1) << "Number of start atom: " << Walker->nr << endl;
+  double Storage[2];
+  atom* Opt_Candidate = NULL;
+  Storage[0]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[1]=999999.;  // This will be an angle looking for the third point.
+cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[Walker->nr] << endl;
 
   for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
     {
-      Walker2 = mol->start;
-      Walker2 = Walker2->next;
-      while (Walker2->nr != (mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr) )
-	{
-	  Walker2 = Walker2->next;
-	}
-
-      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
-    }
-
-  Walker2 = mol->start;
-
-  while (Walker2->nr != int(Storage[0]))
-    {
-      Walker2 = Walker2->next;
-    }
-
+	  Walker2 = mol->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
+    }
+
+  Walker2 = Opt_Candidate;
+  Opt_Candidate=NULL;
   helper.CopyVector(&(Walker->x));
   helper.SubtractVector(&(Walker2->x));
   Oben.ProjectOntoPlane(&helper);
   helper.VectorProduct(&Oben);
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[2]= -10.;  // This will be an angle looking for the third point.
+  Storage[0]=-2.;   	// This will indicate the quadrant.
+  Storage[1]= 9999999.; // This will be an angle looking for the third point.
 
   Chord.CopyVector(&(Walker->x)); // bring into calling function
   Chord.SubtractVector(&(Walker2->x));
 
-printf("Looking for third point candidates \n");
-       Find_next_suitable_point(Walker, Walker2, (mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][0]->rightatom : mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom), 0, &Chord, &helper, &Oben, Storage, RADIUS, mol);
-  Walker3 = mol->start;
-  while (Walker3->nr != int(Storage[0]))
-    {
-      Walker3 = Walker3->next;
-    }
-
-  //Starting Triangle is Walker, Walker2, index Storage[0]
+  cout << Verbose(1) << "Looking for third point candidates \n";
+  Find_next_suitable_point(Walker, Walker2, mol->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
+
+
+  //Starting Triangle is Walker, Walker2, Opt_Candidate
 
   AddPoint(Walker);
   AddPoint(Walker2);
-  AddPoint(Walker3);
+  AddPoint(Opt_Candidate);
+
+  cout << Verbose(1) << "The found starting triangle consists of " << *Walker << ", " << *Walker2 << " and " << *Opt_Candidate << "." << endl;
 
   BPS[0] = new class BoundaryPointSet(Walker);
@@ -1308,5 +1329,5 @@
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
-  BPS[1] = new class BoundaryPointSet(Walker3);
+  BPS[1] = new class BoundaryPointSet(Opt_Candidate);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
@@ -1314,5 +1335,5 @@
   BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  Tesselation::BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
   TrianglesOnBoundaryCount++;
@@ -1324,27 +1345,35 @@
     };
 
-       BTS->GetNormalVector(*BTS->NormalVector);
-
-       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
+       BTS->GetNormalVector(BTS->NormalVector);
+
+       if( BTS->NormalVector.ScalarProduct(&Oben)<0)
 	 {
-	   BTS->NormalVector->Scale(-1);
+	   BTS->NormalVector.Scale(-1);
 	 }
 };
 
 
-void Find_non_convex_border(Tesselation Tess, molecule* mol)
-{
-	printf("Entering finding of non convex hull. \n");
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol)
+{
+	cout << Verbose(1) << "Entering finding of non convex hull. " << endl;
+	cout << flush;
   const double RADIUS =6;
+  LineMap::iterator baseline;
   Tess.Find_starting_triangle(mol, RADIUS);
 
-  for (LineMap::iterator baseline = Tess.LinesOnBoundary.begin(); baseline != Tess.LinesOnBoundary.end(); baseline++)
+  baseline = Tess.LinesOnBoundary.begin();
+  while (baseline != Tess.LinesOnBoundary.end()) {
     if (baseline->second->TrianglesCount == 1)
       {
-		Tess.Find_next_suitable_triangle(mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+    	cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+		Tess.Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+    	cout << Verbose(1) << "End of Tesselation ... " << endl;
       }
     else
       {
-	printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+	cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
       }
-};
+  baseline++;
+  }
+
+};
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision f683feb0f2b4232b7bc42f3893078413c798f659)
+++ src/boundary.hpp	(revision a8bcea6105d46c2f1e551ab69003f9d1cc03e44f)
@@ -36,7 +36,7 @@
     BoundaryPointSet(atom *Walker);
     ~BoundaryPointSet();
-    
+
     void AddLine(class BoundaryLineSet *line);
-    
+
     LineMap lines;
     int LinesCount;
@@ -64,10 +64,10 @@
     BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
     ~BoundaryTriangleSet();
-    
+
     void GetNormalVector(Vector &NormalVector);
-    
+
     class BoundaryPointSet *endpoints[3];
     class BoundaryLineSet *lines[3];
-    Vector *NormalVector;
+    Vector NormalVector;
     int Nr;
 };
@@ -75,14 +75,14 @@
 class Tesselation {
   public:
-    
+
     Tesselation();
     ~Tesselation();
-    
+
     void TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol);
     void GuessStartingTriangle(ofstream *out);
     void AddPoint(atom * Walker);
     void Find_starting_triangle(molecule* mol, const double RADIUS);
-    void Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS);
-    
+    void Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS);
+
     PointMap PointsOnBoundary;
     LineMap LinesOnBoundary;
@@ -105,6 +105,6 @@
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule *mol);
-void Find_non_convex_border(Tesselation Tess, molecule* mol);
+void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule *mol, bool problem);
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol);
 
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision f683feb0f2b4232b7bc42f3893078413c798f659)
+++ src/builder.cpp	(revision a8bcea6105d46c2f1e551ab69003f9d1cc03e44f)
@@ -1,4 +1,4 @@
 /** \file builder.cpp
- * 
+ *
  * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
  * The output is the complete configuration file for PCP for direct use.
@@ -6,43 +6,43 @@
  * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
  * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
- *    
+ *
  */
 
 /*! \mainpage Molecuilder - a molecular set builder
- * 
+ *
  * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
- * 
+ *
  * \section about About the Program
- * 
+ *
  *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
  *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
  *  already constructed atoms.
- * 
+ *
  *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
  *  molecular dynamics implementation.
- * 
+ *
  * \section install Installation
- * 
+ *
  *  Installation should without problems succeed as follows:
  *  -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
  *  -# make
  *  -# make install
- * 
+ *
  *  Further useful commands are
  *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
- *  -# make doxygen-doc: Creates these html pages out of the documented source 
- * 
+ *  -# make doxygen-doc: Creates these html pages out of the documented source
+ *
  * \section run Running
- * 
+ *
  *  The program can be executed by running: ./molecuilder
- * 
+ *
  *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
  *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
- *  later re-execution. 
- * 
+ *  later re-execution.
+ *
  * \section ref References
- *  
+ *
  *  For the special configuration file format, see the documentation of pcp.
- * 
+ *
  */
 
@@ -80,5 +80,5 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
       case 'a': // absolute coordinates of atom
@@ -89,5 +89,5 @@
         mol->AddAtom(first);  // add to molecule
         break;
-        
+
       case 'b': // relative coordinates of atom wrt to reference point
         first = new atom;
@@ -105,5 +105,5 @@
         mol->AddAtom(first);  // add to molecule
         break;
-        
+
       case 'c': // relative coordinates of atom wrt to already placed atom
         first = new atom;
@@ -111,5 +111,5 @@
         do {
           if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
-          second = mol->AskAtom("Enter atom number: ");                
+          second = mol->AskAtom("Enter atom number: ");
           cout << Verbose(0) << "Enter relative coordinates." << endl;
           first->x.AskPosition(mol->cell_size, false);
@@ -121,5 +121,5 @@
         mol->AddAtom(first);  // add to molecule
         break;
-      
+
       case 'd': // two atoms, two angles and a distance
         first = new atom;
@@ -152,5 +152,5 @@
           x.Copyvector(&fourth->x);
           x.SubtractVector(&third->x);
-          
+
           if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
             cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
@@ -167,14 +167,14 @@
           cout << "x: ",
           x.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           z.MakeNormalVector(&second->x,&third->x,&fourth->x);
           cout << "z: ",
           z.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           y.MakeNormalVector(&x,&z);
           cout << "y: ",
           y.Output((ofstream *)&cout);
-          cout << endl;          
-         
+          cout << endl;
+
           // rotate vector around first angle
           first->x.CopyVector(&x);
@@ -182,5 +182,5 @@
           cout << "Rotated vector: ",
           first->x.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           // remove the projection onto the rotation plane of the second angle
           n.CopyVector(&y);
@@ -188,19 +188,19 @@
           cout << "N1: ",
           n.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           first->x.SubtractVector(&n);
           cout << "Subtracted vector: ",
           first->x.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           n.CopyVector(&z);
           n.Scale(first->x.Projection(&z));
           cout << "N2: ",
           n.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           first->x.SubtractVector(&n);
           cout << "2nd subtracted vector: ",
           first->x.Output((ofstream *)&cout);
-          cout << endl;          
-          
+          cout << endl;
+
           // rotate another vector around second angle
           n.CopyVector(&y);
@@ -208,12 +208,12 @@
           cout << "2nd Rotated vector: ",
           n.Output((ofstream *)&cout);
-          cout << endl;          
-          
+          cout << endl;
+
           // add the two linear independent vectors
           first->x.AddVector(&n);
-          first->x.Normalize();         
+          first->x.Normalize();
           first->x.Scale(a);
           first->x.AddVector(&second->x);
-          
+
           cout << Verbose(0) << "resulting coordinates: ";
           first->x.Output((ofstream *)&cout);
@@ -241,5 +241,5 @@
         } while ((j != -1) && (i<128));
         if (i >= 2) {
-          first->x.LSQdistance(atoms, i);             
+          first->x.LSQdistance(atoms, i);
 
           first->x.Output((ofstream *)&cout);
@@ -261,5 +261,5 @@
   Vector x, y;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
   cout << Verbose(0) << " a - on origin" << endl;
@@ -271,5 +271,5 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -327,5 +327,5 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -346,8 +346,8 @@
       second = mol->AskAtom("Enter second atom: ");
 
-      n.CopyVector((const Vector *)&first->x);  
-      n.SubtractVector((const Vector *)&second->x);  
+      n.CopyVector((const Vector *)&first->x);
+      n.SubtractVector((const Vector *)&second->x);
       n.Normalize();
-      break;        
+      break;
     case 'd':
       char shorthand[4];
@@ -363,5 +363,5 @@
         x.x[i] = gsl_vector_get(param.x,i);
         n.x[i] = gsl_vector_get(param.x,i+NDIM);
-      } 
+      }
       gsl_vector_free(param.x);
       cout << Verbose(0) << "Offset vector: ";
@@ -369,5 +369,5 @@
       cout << Verbose(0) << endl;
       n.Normalize();
-      break;       
+      break;
   };
   cout << Verbose(0) << "Alignment vector: ";
@@ -385,5 +385,5 @@
   Vector n;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
   cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
@@ -394,5 +394,5 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -413,8 +413,8 @@
       second = mol->AskAtom("Enter second atom: ");
 
-      n.CopyVector((const Vector *)&first->x);  
-      n.SubtractVector((const Vector *)&second->x);  
+      n.CopyVector((const Vector *)&first->x);
+      n.SubtractVector((const Vector *)&second->x);
       n.Normalize();
-      break;          
+      break;
   };
   cout << Verbose(0) << "Normal vector: ";
@@ -433,5 +433,5 @@
   double tmp1, tmp2;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
   cout << Verbose(0) << " a - state atom for removal by number" << endl;
@@ -442,5 +442,5 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -475,5 +475,5 @@
           mol->RemoveAtom(first);
       }
-      break;          
+      break;
   };
   //mol->Output((ofstream *)&cout);
@@ -492,5 +492,5 @@
   int Z;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
   cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
@@ -514,6 +514,6 @@
       for (int i=MAX_ELEMENTS;i--;)
         min[i] = 0.;
-        
-      second = mol->start;   
+
+      second = mol->start;
       while ((second->next != mol->end)) {
         second = second->next; // advance
@@ -526,10 +526,10 @@
         }
         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
-        //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;         
+        //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
       }
       for (int i=MAX_ELEMENTS;i--;)
         if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
       break;
-      
+
     case 'b':
       first = mol->AskAtom("Enter first atom: ");
@@ -556,5 +556,5 @@
       y.SubtractVector((const Vector *)&second->x);
       cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
-      cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;         
+      cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
       break;
     case 'd':
@@ -600,5 +600,5 @@
   int Order1;
   clock_t start, end;
-  
+
   cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
   cout << Verbose(0) << "What's the desired bond order: ";
@@ -609,5 +609,5 @@
     end = clock();
     cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
-  } else 
+  } else
     cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
 };
@@ -623,5 +623,5 @@
   atom *Walker = mol->start;
   int i, comp, counter=0;
-  
+
   // generate some KeySets
   cout << "Generating KeySets." << endl;
@@ -637,5 +637,5 @@
   cout << "Testing insertion of already present item in KeySets." << endl;
   KeySetTestPair test;
-  test = TestSets[mol->AtomCount-1].insert(Walker->nr); 
+  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
   if (test.second) {
     cout << Verbose(1) << "Insertion worked?!" << endl;
@@ -646,5 +646,5 @@
   TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
 
-  // constructing Graph structure  
+  // constructing Graph structure
   cout << "Generating Subgraph class." << endl;
   Graph Subgraphs;
@@ -657,5 +657,5 @@
   cout << "Testing insertion of already present item in Subgraph." << endl;
   GraphTestPair test2;
-  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 
+  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
   if (test2.second) {
     cout << Verbose(1) << "Insertion worked?!" << endl;
@@ -663,5 +663,5 @@
     cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
   }
-  
+
   // show graphs
   cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
@@ -674,5 +674,5 @@
       if ((*key) > comp)
         cout << (*key) << " ";
-      else 
+      else
         cout << (*key) << "! ";
       comp = (*key);
@@ -716,5 +716,5 @@
   else
     cout << "failed." << endl;
-  
+
   // and save to xyz file
   if (ConfigFileName != NULL) {
@@ -727,5 +727,5 @@
     strcat(filename, ".xyz");
     output.open(filename, ios::trunc);
-  } 
+  }
   cout << Verbose(0) << "Saving of XYZ file ";
   if (mol->MDSteps <= 1) {
@@ -742,5 +742,5 @@
   output.close();
   output.clear();
-  
+
   // and save as MPQC configuration
   if (ConfigFileName != NULL)
@@ -753,5 +753,5 @@
   else
     cout << "failed." << endl;
-  
+
   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
     cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
@@ -785,5 +785,5 @@
   int argptr;
   PathToDatabases = LocalPath;
-  
+
   if (argc > 1) { // config file specified as option
     // 1. : Parse options that just set variables or print help
@@ -798,7 +798,8 @@
           case '?':
             cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
-            cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 
+            cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
             cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
             cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
+            cout << "\t-A <source>\tCreate adjacenzy list from bonds parsed from 'dbond'-style file." <<endl;
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
@@ -812,12 +813,12 @@
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
             cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
-            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
+            cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
-            cout << "\t-T <file> Store temperatures from the config file in <file>." << endl; 
+            cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
-            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 
+            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
             cout << "\t-v/-V\t\tGives version information." << endl;
             cout << "Note: config files must not begin with '-' !" << endl;
@@ -854,5 +855,5 @@
         argptr++;
     } while (argptr < argc);
-    
+
     // 2. Parse the element database
     if (periode->LoadPeriodentafel(PathToDatabases)) {
@@ -863,5 +864,5 @@
       return 1;
     }
-    
+
     // 3. Find config file name and parse if possible
     if (argv[1][0] != '-') {
@@ -902,5 +903,5 @@
     } else
       config_present = absent;
-    
+
     // 4. parse again through options, now for those depending on elements db and config presence
     argptr = 1;
@@ -946,5 +947,5 @@
                     config_present = present;
                 } else
-                  cerr << Verbose(1) << "Could not find the specified element." << endl; 
+                  cerr << Verbose(1) << "Could not find the specified element." << endl;
                 argptr+=4;
               }
@@ -956,5 +957,5 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
-            case 'D': 
+            case 'D':
               ExitFlag = 1;
               {
@@ -1002,4 +1003,17 @@
               }
               break;
+            case 'A':
+            	ExitFlag = 1;
+            	if ((argptr >= argc) || (argv[argptr][0] == '-')){
+            		ExitFlag =255;
+            		cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
+            	}
+            	else{
+            		cout << "Parsing bonds from " << argv[argptr] << "." << endl;
+                    ifstream *input = new ifstream(argv[argptr]);
+            		mol->CreateAdjacencyList2((ofstream *)&cout, input);
+            		input->close();
+            	}
+            	break;
             case 'T':
               ExitFlag = 1;
@@ -1168,11 +1182,12 @@
             case 'o':
               ExitFlag = 1;
-              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+              if ((argptr >= argc) || (argv[argptr][0] == '-')){
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
               } else {
                 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+                ofstream *output = new ofstream(argv[argptr], ios::trunc);
+                //$$$
                 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
-                ofstream *output = new ofstream(argv[argptr], ios::trunc);
                 VolumeOfConvexEnvelope((ofstream *)&cout, output, &configuration, NULL, mol);
                 output->close();
@@ -1269,5 +1284,5 @@
                       mol->Translate(&x);
                     }
-                    mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+                    mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
                   }
                 }
@@ -1332,5 +1347,5 @@
   if (j == 1) return 0; // just for -v and -h options
   if (j) return j;  // something went wrong
-  
+
   // General stuff
   if (mol->cell_size[0] == 0.) {
@@ -1346,5 +1361,5 @@
   // now the main construction loop
   cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
-  do {    
+  do {
     cout << Verbose(0) << endl << endl;
     cout << Verbose(0) << "============Element list=======================" << endl;
@@ -1365,6 +1380,6 @@
     cout << Verbose(0) << "-----------------------------------------------" << endl;
     cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
-    cout << Verbose(0) << "i - realign molecule" << endl; 
-    cout << Verbose(0) << "m - mirror all molecules" << endl; 
+    cout << Verbose(0) << "i - realign molecule" << endl;
+    cout << Verbose(0) << "m - mirror all molecules" << endl;
     cout << Verbose(0) << "t - translate molecule by vector" << endl;
     cout << Verbose(0) << "c - scale by unit transformation" << endl;
@@ -1377,12 +1392,12 @@
     cout << Verbose(0) << "Input: ";
     cin >> choice;
-    
+
     switch (choice) {
       default:
       case 'a': // add atom
         AddAtoms(periode, mol);
-        choice = 'a'; 
-        break;
-      
+        choice = 'a';
+        break;
+
       case 'b': // scale a bond
         cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
@@ -1400,8 +1415,8 @@
         }
         //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
-        //second->Output(second->type->No, 1, (ofstream *)&cout);       
-        break;
-
-      case 'c': // unit scaling of the metric 
+        //second->Output(second->type->No, 1, (ofstream *)&cout);
+        break;
+
+      case 'c': // unit scaling of the metric
        cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
        cout << Verbose(0) << "Enter three factors: ";
@@ -1414,5 +1429,5 @@
        delete[](factor);
        break;
-       
+
       case 'd': // duplicate the periodic cell along a given axis, given times
         cout << Verbose(0) << "State the axis [(+-)123]: ";
@@ -1420,5 +1435,5 @@
         cout << Verbose(0) << "State the factor: ";
         cin >> faktor;
-        
+
         mol->CountAtoms((ofstream *)&cout);  // recount atoms
         if (mol->AtomCount != 0) {  // if there is more than none
@@ -1461,21 +1476,21 @@
             mol->Translate(&x);
           }
-          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
         }
         break;
-      
+
       case 'e': // edit each field of the configuration
        configuration.Edit(mol);
        break;
-  
+
       case 'f':
         FragmentAtoms(mol, &configuration);
         break;
-        
+
       case 'g': // center the atoms
         CenterAtoms(mol);
         break;
-        
-      case 'i': // align all atoms 
+
+      case 'i': // align all atoms
         AlignAtoms(periode, mol);
         break;
@@ -1488,5 +1503,5 @@
         MirrorAtoms(mol);
         break;
-       
+
       case 'o': // create the connection matrix
         {
@@ -1509,5 +1524,5 @@
         }
         break;
-        
+
       case 'p': // parse and XYZ file
         cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
@@ -1518,11 +1533,11 @@
         break;
 
-      case 'q': // quit 
-        break;
-        
+      case 'q': // quit
+        break;
+
       case 'r': // remove atom
-        RemoveAtoms(mol);        
-        break;
-        
+        RemoveAtoms(mol);
+        break;
+
       case 's': // save to config file
         SaveConfig(ConfigFileName, &configuration, periode, mol);
@@ -1530,13 +1545,13 @@
 
       case 't': // translate all atoms
-       cout << Verbose(0) << "Enter translation vector." << endl;       
+       cout << Verbose(0) << "Enter translation vector." << endl;
        x.AskPosition(mol->cell_size,0);
        mol->Translate((const Vector *)&x);
        break;
-  
+
       case 'T':
         testroutine(mol);
         break;
-      
+
       case 'u': // change an atom's element
         first = NULL;
@@ -1545,12 +1560,12 @@
           cin >> Z;
         } while ((first = mol->FindAtom(Z)) == NULL);
-        cout << Verbose(0) << "New element by atomic number Z: ";        
+        cout << Verbose(0) << "New element by atomic number Z: ";
         cin >> Z;
         first->type = periode->FindElement(Z);
-        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;    
+        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
         break;
     };
   } while (choice != 'q');
-  
+
   // save element data base
   if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
