Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/atom.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -13,5 +13,5 @@
 atom::atom() 
 {
-	Name = NULL;
+  Name = NULL;
   previous = NULL;
   next = NULL;
@@ -35,5 +35,5 @@
 atom::~atom() 
 {
-	Free((void **)&Name, "atom::~atom: *Name");
+  Free((void **)&Name, "atom::~atom: *Name");
   Free((void **)&ComponentNr, "atom::~atom: *ComponentNr");
 };
Index: src/bond.cpp
===================================================================
--- src/bond.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/bond.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -14,11 +14,11 @@
 bond::bond()
 {
-	leftatom = NULL;
-	rightatom = NULL;
+  leftatom = NULL;
+  rightatom = NULL;
   previous = NULL;
   next = NULL;
-	nr = -1;
-	HydrogenBond = 0;
-	BondDegree = 0;
+  nr = -1;
+  HydrogenBond = 0;
+  BondDegree = 0;
   Used = white;
   Cyclic = false;
@@ -34,14 +34,14 @@
 bond::bond(atom *left, atom *right, int degree=1, int number=0)
 {
-	leftatom = left;
-	rightatom = right;
+  leftatom = left;
+  rightatom = right;
   previous = NULL;
   next = NULL;
-	HydrogenBond = 0;
+  HydrogenBond = 0;
   if ((left != NULL) && (right != NULL)) {
-  	if ((left->type != NULL) && (left->type->Z == 1))
-  		HydrogenBond++;
-  	if ((right->type != NULL) && (right->type->Z == 1))
-  		HydrogenBond++;
+    if ((left->type != NULL) && (left->type->Z == 1))
+      HydrogenBond++;
+    if ((right->type != NULL) && (right->type->Z == 1))
+      HydrogenBond++;
   }
   BondDegree = degree;
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/boundary.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -86,19 +86,19 @@
 BoundaryLineSet::~BoundaryLineSet()
 {
-	for (int i = 0; i < 2; i++) {
-		cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-		endpoints[i]->lines.erase(Nr);
-		LineMap::iterator tester = endpoints[i]->lines.begin();
-		tester++;
-		if (tester == endpoints[i]->lines.end()) {
-			cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-			if (endpoints[i] != NULL) {
-				delete(endpoints[i]);
-				endpoints[i] = NULL;
-			} else
-				cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
-		} else
-			cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
-	}
+  for (int i = 0; i < 2; i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+    endpoints[i]->lines.erase(Nr);
+    LineMap::iterator tester = endpoints[i]->lines.begin();
+    tester++;
+    if (tester == endpoints[i]->lines.end()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      if (endpoints[i] != NULL) {
+        delete(endpoints[i]);
+        endpoints[i] = NULL;
+      } else
+        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
 }
 ;
@@ -181,17 +181,17 @@
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
-	for (int i = 0; i < 3; i++) {
-		cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
-		lines[i]->triangles.erase(Nr);
-		if (lines[i]->triangles.empty()) {
-			cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-			if (lines[i] != NULL) {
-				delete (lines[i]);
-				lines[i] = NULL;
-			} else
-				cerr << "ERROR: This line " << i << " has already been free'd." << endl;
-		} else
-			cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
-	}
+  for (int i = 0; i < 3; i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    if (lines[i]->triangles.empty()) {
+      cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+      if (lines[i] != NULL) {
+        delete (lines[i]);
+        lines[i] = NULL;
+      } else
+        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+  }
 }
 ;
@@ -573,48 +573,48 @@
 void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
 {
-	atom *Walker = mol->start;
-	bond *Binder = mol->first;
-	int i;
-	Vector *center = mol->DetermineCenterOfAll(out);
-	if (rasterfile != NULL) {
-		//cout << Verbose(1) << "Writing Raster3D file ... ";
-		*rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
-		*rasterfile << "@header.r3d" << endl;
-		*rasterfile << "# All atoms as spheres" << endl;
-		while (Walker->next != mol->end) {
-			Walker = Walker->next;
-			*rasterfile << "2" << endl << "  ";	// 2 is sphere type
-			for (i=0;i<NDIM;i++)
-				*rasterfile << Walker->x.x[i]+center->x[i] << " ";
-			*rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
-		}
-
-		*rasterfile << "# All bonds as vertices" << endl;
-		while (Binder->next != mol->last) {
-			Binder = Binder->next;
-			*rasterfile << "3" << endl << "  ";	// 2 is round-ended cylinder type
-			for (i=0;i<NDIM;i++)
-				*rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
-			*rasterfile << "\t0.03\t";
-			for (i=0;i<NDIM;i++)
-				*rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
-			*rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
-		}
-
-		*rasterfile << "# All tesselation triangles" << endl;
-		for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-			*rasterfile << "1" << endl << "  ";	// 1 is triangle type
-			for (i=0;i<3;i++) {	// print each node
-				for (int j=0;j<NDIM;j++)	// and for each node all NDIM coordinates
-					*rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
-				*rasterfile << "\t";
-			}
-			*rasterfile << "1. 0. 0." << endl;	// red as colour
-			*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;	// 18 is transparency type for previous object
-		}
-	} else {
-		cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
-	}
-	delete(center);
+  atom *Walker = mol->start;
+  bond *Binder = mol->first;
+  int i;
+  Vector *center = mol->DetermineCenterOfAll(out);
+  if (rasterfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+    *rasterfile << "@header.r3d" << endl;
+    *rasterfile << "# All atoms as spheres" << endl;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Walker->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
+
+    *rasterfile << "# All bonds as vertices" << endl;
+    while (Binder->next != mol->last) {
+      Binder = Binder->next;
+      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t";
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
+
+    *rasterfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
+      for (i=0;i<3;i++) {  // print each node
+        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+        *rasterfile << "\t";
+      }
+      *rasterfile << "1. 0. 0." << endl;  // red as colour
+      *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  delete(center);
 };
 
@@ -946,26 +946,26 @@
 Tesselation::~Tesselation()
 {
-	cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-	for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
-		if (runner->second != NULL) {
-			delete (runner->second);
-			runner->second = NULL;
-		} else
-			cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
-	}
-	for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {
-		if (runner->second != NULL) {
-			delete (runner->second);
-			runner->second = NULL;
-		} else
-			cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;
-	}
-	for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-		if (runner->second != NULL) {
-			delete (runner->second);
-			runner->second = NULL;
-		} else
-			cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;
-	}
+  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+  }
+  for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;
+  }
+  for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;
+  }
 }
 ;
@@ -1585,156 +1585,156 @@
     atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
 {
-	cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-	cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
-	cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
-	cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
-	cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
-	cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
-	/* OldNormal is normal vector on the old triangle
-	 * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
-	 * Chord points from b to a!!!
-	 */
-	Vector dif_a; //Vector from a to candidate
-	Vector dif_b; //Vector from b to candidate
-	Vector AngleCheck;
-	Vector TempNormal, Umkreismittelpunkt;
-	Vector Mittelpunkt;
-
-	double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
-	double BallAngle, AlternativeSign;
-	atom *Walker; // variable atom point
-
-	Vector NewUmkreismittelpunkt;
-
-	if (a != Candidate and b != Candidate and c != Candidate) {
-		cout << Verbose(3) << "We have a unique candidate!" << endl;
-		dif_a.CopyVector(&(a->x));
-		dif_a.SubtractVector(&(Candidate->x));
-		dif_b.CopyVector(&(b->x));
-		dif_b.SubtractVector(&(Candidate->x));
-		AngleCheck.CopyVector(&(Candidate->x));
-		AngleCheck.SubtractVector(&(a->x));
-		AngleCheck.ProjectOntoPlane(Chord);
-
-		SideA = dif_b.Norm();
-		SideB = dif_a.Norm();
-		SideC = Chord->Norm();
-		//Chord->Scale(-1);
-
-		alpha = Chord->Angle(&dif_a);
-		beta = M_PI - Chord->Angle(&dif_b);
-		gamma = dif_a.Angle(&dif_b);
-
-		cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-
-		if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
-			cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-			cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-		}
-
-		if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
-			Umkreisradius = SideA / 2.0 / sin(alpha);
-			//cout << Umkreisradius << endl;
-			//cout << SideB / 2.0 / sin(beta) << endl;
-			//cout << SideC / 2.0 / sin(gamma) << endl;
-
-			if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
-				cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
-				cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
-				sign = AngleCheck.ScalarProduct(direction1);
-				if (fabs(sign)<MYEPSILON) {
-					if (AngleCheck.ScalarProduct(OldNormal)<0) {
-						sign =0;
-						AlternativeSign=1;
-					} else {
-						sign =0;
-						AlternativeSign=-1;
-					}
-				} else {
-					sign /= fabs(sign);
-					cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
-				}
-
-				Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
-
-				AngleCheck.CopyVector(&ReferencePoint);
-				AngleCheck.Scale(-1);
-				//cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-				AngleCheck.AddVector(&Mittelpunkt);
-				//cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-				cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
-
-				BallAngle = AngleCheck.Angle(OldNormal);
-				cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-
-				//cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
-				//cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-
-				cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
-
-				NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
-
-				if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
-					if (Storage[0]< -1.5) { // first Candidate at all
-						if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
-							cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
-							Opt_Candidate = Candidate;
-							Storage[0] = sign;
-							Storage[1] = AlternativeSign;
-							Storage[2] = BallAngle;
-							cout << " angle " << Storage[2] << " and Up/Down "
-							<< Storage[0] << endl;
-						} else
-							cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
-					} else {
-						if ( Storage[2] > BallAngle) {
-							if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
-								cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
-								Opt_Candidate = Candidate;
-								Storage[0] = sign;
-								Storage[1] = AlternativeSign;
-								Storage[2] = BallAngle;
-								cout << " angle " << Storage[2] << " and Up/Down "
-								<< Storage[0] << endl;
-							} else
-								cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
-						} else {
-							if (DEBUG) {
-								cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
-							}
-						}
-					}
-				} else {
-					if (DEBUG) {
-						cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
-					}
-				}
-			} else {
-				if (DEBUG) {
-					cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
-				}
-			}
-		} else {
-			if (DEBUG) {
-				cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
-			}
-		}
-	} else {
-		if (DEBUG) {
-			cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
-		}
-	}
-
-	if (RecursionLevel < 5) { // Seven is the recursion level threshold.
-		for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
-			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-			if (Walker == Parent) { // don't go back the same bond
-				continue;
-			} else {
-				Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
-			}
-		}
-	}
-	cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
+  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
+  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
+  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
+  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
+  /* OldNormal is normal vector on the old triangle
+   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck;
+  Vector TempNormal, Umkreismittelpunkt;
+  Vector Mittelpunkt;
+
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
+  double BallAngle, AlternativeSign;
+  atom *Walker; // variable atom point
+
+  Vector NewUmkreismittelpunkt;
+
+  if (a != Candidate and b != Candidate and c != Candidate) {
+    cout << Verbose(3) << "We have a unique candidate!" << endl;
+    dif_a.CopyVector(&(a->x));
+    dif_a.SubtractVector(&(Candidate->x));
+    dif_b.CopyVector(&(b->x));
+    dif_b.SubtractVector(&(Candidate->x));
+    AngleCheck.CopyVector(&(Candidate->x));
+    AngleCheck.SubtractVector(&(a->x));
+    AngleCheck.ProjectOntoPlane(Chord);
+
+    SideA = dif_b.Norm();
+    SideB = dif_a.Norm();
+    SideC = Chord->Norm();
+    //Chord->Scale(-1);
+
+    alpha = Chord->Angle(&dif_a);
+    beta = M_PI - Chord->Angle(&dif_b);
+    gamma = dif_a.Angle(&dif_b);
+
+    cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+
+    if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
+      cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+      cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+    }
+
+    if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
+      Umkreisradius = SideA / 2.0 / sin(alpha);
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
+
+      if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
+        cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
+        cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
+        sign = AngleCheck.ScalarProduct(direction1);
+        if (fabs(sign)<MYEPSILON) {
+          if (AngleCheck.ScalarProduct(OldNormal)<0) {
+            sign =0;
+            AlternativeSign=1;
+          } else {
+            sign =0;
+            AlternativeSign=-1;
+          }
+        } else {
+          sign /= fabs(sign);
+          cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
+        }
+
+        Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+
+        AngleCheck.CopyVector(&ReferencePoint);
+        AngleCheck.Scale(-1);
+        //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+        AngleCheck.AddVector(&Mittelpunkt);
+        //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+        cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
+
+        BallAngle = AngleCheck.Angle(OldNormal);
+        cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+
+        //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+        //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+        cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+
+        NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
+
+        if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
+          if (Storage[0]< -1.5) { // first Candidate at all
+            if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+              cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
+              Opt_Candidate = Candidate;
+              Storage[0] = sign;
+              Storage[1] = AlternativeSign;
+              Storage[2] = BallAngle;
+              cout << " angle " << Storage[2] << " and Up/Down "
+              << Storage[0] << endl;
+            } else
+              cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+          } else {
+            if ( Storage[2] > BallAngle) {
+              if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+                cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
+                Opt_Candidate = Candidate;
+                Storage[0] = sign;
+                Storage[1] = AlternativeSign;
+                Storage[2] = BallAngle;
+                cout << " angle " << Storage[2] << " and Up/Down "
+                << Storage[0] << endl;
+              } else
+                cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+            } else {
+              if (DEBUG) {
+                cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
+              }
+            }
+          }
+        } else {
+          if (DEBUG) {
+            cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
+          }
+        }
+      } else {
+        if (DEBUG) {
+          cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
+        }
+      }
+    } else {
+      if (DEBUG) {
+        cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
+      }
+    }
+  } else {
+    if (DEBUG) {
+      cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
+    }
+  }
+
+  if (RecursionLevel < 5) { // Seven is the recursion level threshold.
+    for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+      Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+      if (Walker == Parent) { // don't go back the same bond
+        continue;
+      } else {
+        Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
+      }
+    }
+  }
+  cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
 }
 ;
@@ -1767,142 +1767,142 @@
     atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
 {
-	/* OldNormal is normal vector on the old triangle
-	 * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
-	 * Chord points from b to a!!!
-	 */
-	Vector dif_a; //Vector from a to candidate
-	Vector dif_b; //Vector from b to candidate
-	Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
-	Vector TempNormal, Umkreismittelpunkt, Mittelpunkt, helper;
-
-	double CurrentEpsilon = 0.1;
-	double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
-	double BallAngle;
-	atom *Walker; // variable atom point
-
-
-	dif_a.CopyVector(&(a->x));
-	dif_a.SubtractVector(&(Candidate->x));
-	dif_b.CopyVector(&(b->x));
-	dif_b.SubtractVector(&(Candidate->x));
-	DirectionCheckPoint.CopyVector(&dif_a);
-	DirectionCheckPoint.Scale(-1);
-	DirectionCheckPoint.ProjectOntoPlane(Chord);
-
-	SideA = dif_b.Norm();
-	SideB = dif_a.Norm();
-	SideC = Chord->Norm();
-	//Chord->Scale(-1);
-
-	alpha = Chord->Angle(&dif_a);
-	beta = M_PI - Chord->Angle(&dif_b);
-	gamma = dif_a.Angle(&dif_b);
-
-
-	if (DEBUG) {
-		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) {
-		//      alpha = dif_a.Angle(&dif_b) / 2.;
-		//      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
-		//      SideB = dif_a.Norm();
-		//      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
-		//          alpha); // note this is squared of center line length
-		//      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-		// Those are remains from Freddie. Needed?
-
-		Umkreisradius = SideA / 2.0 / sin(alpha);
-		//cout << Umkreisradius << endl;
-		//cout << SideB / 2.0 / sin(beta) << endl;
-		//cout << SideC / 2.0 / sin(gamma) << endl;
-
-		if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
-			// intermediate calculations to aquire centre of sphere, called Mittelpunkt:
-			Umkreismittelpunkt.Zero();
-			helper.CopyVector(&a->x);
-			helper.Scale(sin(2.*alpha));
-			Umkreismittelpunkt.AddVector(&helper);
-			helper.CopyVector(&b->x);
-			helper.Scale(sin(2.*beta));
-			Umkreismittelpunkt.AddVector(&helper);
-			helper.CopyVector(&Candidate->x);
-			helper.Scale(sin(2.*gamma));
-			Umkreismittelpunkt.AddVector(&helper);
-			//Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
-			Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
-			TempNormal.CopyVector(&dif_a);
-			TempNormal.VectorProduct(&dif_b);
-			if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
-				TempNormal.Scale(-1);
-			}
-			TempNormal.Normalize();
-			Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-			TempNormal.Scale(Restradius);
-
-			Mittelpunkt.CopyVector(&Umkreismittelpunkt);
-			Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
-
-			AngleCheck.CopyVector(Chord);
-			AngleCheck.Scale(-0.5);
-			AngleCheck.SubtractVector(&(b->x));
-			AngleCheckReference.CopyVector(&AngleCheck);
-			AngleCheckReference.AddVector(Opt_Mittelpunkt);
-			AngleCheck.AddVector(&Mittelpunkt);
-
-			BallAngle = AngleCheck.Angle(&AngleCheckReference);
-
-			d1->ProjectOntoPlane(&AngleCheckReference);
-			sign = AngleCheck.ScalarProduct(d1);
-			if (fabs(sign) < MYEPSILON)
-				sign = 0;
-			else
-				sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-
-
-			if (Storage[0]< -1.5) { // first Candidate at all
-				cout << "Next better candidate is " << *Candidate << " with ";
-				Opt_Candidate = Candidate;
-				Storage[0] = sign;
-				Storage[1] = BallAngle;
-				Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-				cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-			} else {
-				/*
-				 * removed due to change in criterium, now checking angle of ball to old normal.
-				//We will now check for non interference, that is if the new candidate would have the Opt_Candidate
-				//within the ball.
-
-				Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
-				//cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
-
-
-				if (Distance >RADIUS) { // We have no interference and may now check whether the new point is better.
-				 */
-					//cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
-
-				if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))) { //This will give absolute preference to those in "right-hand" quadrants
-					//(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
-					cout << "Next better candidate is " << *Candidate << " with ";
-					Opt_Candidate = Candidate;
-					Storage[0] = sign;
-					Storage[1] = BallAngle;
-					Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-					cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-				} else {
-					if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0 && Storage[1] > BallAngle) || (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0 && Storage[1] < BallAngle)) {
-						//Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-						cout << "Next better candidate is " << *Candidate << " with ";
-						Opt_Candidate = Candidate;
-						Storage[0] = sign;
-						Storage[1] = BallAngle;
-						Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-						cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-					}
-				}
-			}
+  /* OldNormal is normal vector on the old triangle
+   * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
+  Vector TempNormal, Umkreismittelpunkt, Mittelpunkt, helper;
+
+  double CurrentEpsilon = 0.1;
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
+  double BallAngle;
+  atom *Walker; // variable atom point
+
+
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
+  DirectionCheckPoint.CopyVector(&dif_a);
+  DirectionCheckPoint.Scale(-1);
+  DirectionCheckPoint.ProjectOntoPlane(Chord);
+
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+
+
+  if (DEBUG) {
+    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) {
+    //      alpha = dif_a.Angle(&dif_b) / 2.;
+    //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
+    //      SideB = dif_a.Norm();
+    //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
+    //          alpha); // note this is squared of center line length
+    //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
+    // Those are remains from Freddie. Needed?
+
+    Umkreisradius = SideA / 2.0 / sin(alpha);
+    //cout << Umkreisradius << endl;
+    //cout << SideB / 2.0 / sin(beta) << endl;
+    //cout << SideC / 2.0 / sin(gamma) << endl;
+
+    if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
+      // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
+      Umkreismittelpunkt.Zero();
+      helper.CopyVector(&a->x);
+      helper.Scale(sin(2.*alpha));
+      Umkreismittelpunkt.AddVector(&helper);
+      helper.CopyVector(&b->x);
+      helper.Scale(sin(2.*beta));
+      Umkreismittelpunkt.AddVector(&helper);
+      helper.CopyVector(&Candidate->x);
+      helper.Scale(sin(2.*gamma));
+      Umkreismittelpunkt.AddVector(&helper);
+      //Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
+      Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+
+      TempNormal.CopyVector(&dif_a);
+      TempNormal.VectorProduct(&dif_b);
+      if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
+        TempNormal.Scale(-1);
+      }
+      TempNormal.Normalize();
+      Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+      TempNormal.Scale(Restradius);
+
+      Mittelpunkt.CopyVector(&Umkreismittelpunkt);
+      Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
+
+      AngleCheck.CopyVector(Chord);
+      AngleCheck.Scale(-0.5);
+      AngleCheck.SubtractVector(&(b->x));
+      AngleCheckReference.CopyVector(&AngleCheck);
+      AngleCheckReference.AddVector(Opt_Mittelpunkt);
+      AngleCheck.AddVector(&Mittelpunkt);
+
+      BallAngle = AngleCheck.Angle(&AngleCheckReference);
+
+      d1->ProjectOntoPlane(&AngleCheckReference);
+      sign = AngleCheck.ScalarProduct(d1);
+      if (fabs(sign) < MYEPSILON)
+        sign = 0;
+      else
+        sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+
+
+      if (Storage[0]< -1.5) { // first Candidate at all
+        cout << "Next better candidate is " << *Candidate << " with ";
+        Opt_Candidate = Candidate;
+        Storage[0] = sign;
+        Storage[1] = BallAngle;
+        Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+        cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+      } else {
+        /*
+         * removed due to change in criterium, now checking angle of ball to old normal.
+        //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
+        //within the ball.
+
+        Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
+        //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+
+
+        if (Distance >RADIUS) { // We have no interference and may now check whether the new point is better.
+         */
+          //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
+
+        if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))) { //This will give absolute preference to those in "right-hand" quadrants
+          //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
+          cout << "Next better candidate is " << *Candidate << " with ";
+          Opt_Candidate = Candidate;
+          Storage[0] = sign;
+          Storage[1] = BallAngle;
+          Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+          cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+        } else {
+          if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0 && Storage[1] > BallAngle) || (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0 && Storage[1] < BallAngle)) {
+            //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+            cout << "Next better candidate is " << *Candidate << " with ";
+            Opt_Candidate = Candidate;
+            Storage[0] = sign;
+            Storage[1] = BallAngle;
+            Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+            cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+          }
+        }
+      }
 /*
                * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
@@ -1947,25 +1947,25 @@
                 }
                 */
-		} else {
-			if (DEBUG) {
-				cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
-			}
-		}
-	} else {
-		if (DEBUG) {
-			cout << "identisch mit Ursprungslinie" << endl;
-		}
-	}
-
-	if (RecursionLevel < 9) { // Five is the recursion level threshold.
-		for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
-			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-			if (Walker == Parent) { // don't go back the same bond
-				continue;
-			} else {
-				Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel+1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS, mol); //call function again
-			}
-		}
-	}
+    } else {
+      if (DEBUG) {
+        cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+      }
+    }
+  } else {
+    if (DEBUG) {
+      cout << "identisch mit Ursprungslinie" << endl;
+    }
+  }
+
+  if (RecursionLevel < 9) { // Five is the recursion level threshold.
+    for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+      Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+      if (Walker == Parent) { // don't go back the same bond
+        continue;
+      } else {
+        Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel+1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS, mol); //call function again
+      }
+    }
+  }
 };
 
@@ -1984,166 +1984,166 @@
     const double& RADIUS, int N, const char *tempbasename)
 {
-	cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
-	Vector direction1;
-	Vector helper;
-	Vector Chord;
-	ofstream *tempstream = NULL;
-	char NumberName[255];
-	double tmp;
-	//atom* Walker;
-	atom* OldThirdPoint;
-
-	double Storage[3];
-	Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
-	Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
-	Storage[2] = 9999999.;
-	atom* Opt_Candidate = NULL;
-	Vector Opt_Mittelpunkt;
-
-	cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-
-	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) {
-			OldThirdPoint = T.endpoints[i]->node;
-			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);
-	}
-	cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
-
-	Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
-	Chord.SubtractVector(&(Line.endpoints[1]->node->x));
-	cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
-
-
-	Vector Umkreismittelpunkt, a, b, c;
-	double alpha, beta, gamma;
-	a.CopyVector(&(T.endpoints[0]->node->x));
-	b.CopyVector(&(T.endpoints[1]->node->x));
-	c.CopyVector(&(T.endpoints[2]->node->x));
-	a.SubtractVector(&(T.endpoints[1]->node->x));
-	b.SubtractVector(&(T.endpoints[2]->node->x));
-	c.SubtractVector(&(T.endpoints[0]->node->x));
-
-	alpha = M_PI - a.Angle(&c);
-	beta = M_PI - b.Angle(&a);
-	gamma = M_PI - c.Angle(&b);
-	if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-		cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-
-	Umkreismittelpunkt.Zero();
-	helper.CopyVector(&T.endpoints[0]->node->x);
-	helper.Scale(sin(2.*alpha));
-	Umkreismittelpunkt.AddVector(&helper);
-	helper.CopyVector(&T.endpoints[1]->node->x);
-	helper.Scale(sin(2.*beta));
-	Umkreismittelpunkt.AddVector(&helper);
-	helper.CopyVector(&T.endpoints[2]->node->x);
-	helper.Scale(sin(2.*gamma));
-	Umkreismittelpunkt.AddVector(&helper);
-	//Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
-	//cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-	Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-	//cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-	cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
-	cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
-
-	cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-	cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
-	if (DEBUG)
-		cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
-	tmp = 0;
-	for (int i=0;i<NDIM;i++) {
-		helper.CopyVector(&T.endpoints[i]->node->x);
-		helper.SubtractVector(&Umkreismittelpunkt);
-		if (DEBUG)
-			cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
-		if (tmp == 0) // set on first time for comparison against next ones
-			tmp = helper.Norm();
-		if (fabs(helper.Norm() - tmp) > MYEPSILON)
-			cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
-	}
-
-	cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
-
-	Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-	Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-	if (Opt_Candidate == NULL) {
-		cerr << "WARNING: Could not find a suitable candidate." << endl;
-		return false;
-	}
-	cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << endl;
-
-	// check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-	bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node);
-
-	if (flag) { // if so, add
-		AddTrianglePoint(Opt_Candidate, 0);
-		AddTrianglePoint(Line.endpoints[0]->node, 1);
-		AddTrianglePoint(Line.endpoints[1]->node, 2);
-
-		AddTriangleLine(TPS[0], TPS[1], 0);
-		AddTriangleLine(TPS[0], TPS[2], 1);
-		AddTriangleLine(TPS[1], TPS[2], 2);
-
-		BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-		AddTriangleToLines();
-
-		BTS->GetNormalVector(BTS->NormalVector);
-
-		if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  || (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) || (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) ) {
-			BTS->NormalVector.Scale(-1);
-		};
-		cout << Verbose(1) << "New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-		cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
-	} else { // else, yell and do nothing
-		cout << Verbose(1) << "This triangle consisting of ";
-		cout << *Opt_Candidate << ", ";
-		cout << *Line.endpoints[0]->node << " and ";
-		cout << *Line.endpoints[1]->node << " ";
-		cout << "is invalid!" << endl;
-		return false;
-	}
-
-	if ((TrianglesOnBoundaryCount % 10) == 0) {
-		sprintf(NumberName, "-%d", TriangleFilesWritten);
-		if (DoTecplotOutput) {
-			string NameofTempFile(tempbasename);
-			NameofTempFile.append(NumberName);
-			NameofTempFile.append(TecplotSuffix);
-			cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-			tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-			write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
-			tempstream->close();
-			tempstream->flush();
-			delete(tempstream);
-		}
-
-		if (DoRaster3DOutput) {
-			string NameofTempFile(tempbasename);
-			NameofTempFile.append(NumberName);
-			NameofTempFile.append(Raster3DSuffix);
-			cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-			tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-			write_raster3d_file(out, tempstream, this, mol);
-			tempstream->close();
-			tempstream->flush();
-			delete(tempstream);
-		}
-		if (DoTecplotOutput || DoRaster3DOutput)
-			TriangleFilesWritten++;
-	}
-
-	cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
-	return true;
+  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+  Vector direction1;
+  Vector helper;
+  Vector Chord;
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  double tmp;
+  //atom* Walker;
+  atom* OldThirdPoint;
+
+  double Storage[3];
+  Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[2] = 9999999.;
+  atom* Opt_Candidate = NULL;
+  Vector Opt_Mittelpunkt;
+
+  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+
+  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) {
+      OldThirdPoint = T.endpoints[i]->node;
+      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);
+  }
+  cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
+
+  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
+  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
+  cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
+
+
+  Vector Umkreismittelpunkt, a, b, c;
+  double alpha, beta, gamma;
+  a.CopyVector(&(T.endpoints[0]->node->x));
+  b.CopyVector(&(T.endpoints[1]->node->x));
+  c.CopyVector(&(T.endpoints[2]->node->x));
+  a.SubtractVector(&(T.endpoints[1]->node->x));
+  b.SubtractVector(&(T.endpoints[2]->node->x));
+  c.SubtractVector(&(T.endpoints[0]->node->x));
+
+  alpha = M_PI - a.Angle(&c);
+  beta = M_PI - b.Angle(&a);
+  gamma = M_PI - c.Angle(&b);
+  if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
+    cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+
+  Umkreismittelpunkt.Zero();
+  helper.CopyVector(&T.endpoints[0]->node->x);
+  helper.Scale(sin(2.*alpha));
+  Umkreismittelpunkt.AddVector(&helper);
+  helper.CopyVector(&T.endpoints[1]->node->x);
+  helper.Scale(sin(2.*beta));
+  Umkreismittelpunkt.AddVector(&helper);
+  helper.CopyVector(&T.endpoints[2]->node->x);
+  helper.Scale(sin(2.*gamma));
+  Umkreismittelpunkt.AddVector(&helper);
+  //Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+  //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+  Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+  //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+  cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
+  cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
+
+  cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+  cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
+  if (DEBUG)
+    cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
+  tmp = 0;
+  for (int i=0;i<NDIM;i++) {
+    helper.CopyVector(&T.endpoints[i]->node->x);
+    helper.SubtractVector(&Umkreismittelpunkt);
+    if (DEBUG)
+      cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
+    if (tmp == 0) // set on first time for comparison against next ones
+      tmp = helper.Norm();
+    if (fabs(helper.Norm() - tmp) > MYEPSILON)
+      cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
+  }
+
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+  if (Opt_Candidate == NULL) {
+    cerr << "WARNING: Could not find a suitable candidate." << endl;
+    return false;
+  }
+  cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << endl;
+
+  // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+  bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node);
+
+  if (flag) { // if so, add
+    AddTrianglePoint(Opt_Candidate, 0);
+    AddTrianglePoint(Line.endpoints[0]->node, 1);
+    AddTrianglePoint(Line.endpoints[1]->node, 2);
+
+    AddTriangleLine(TPS[0], TPS[1], 0);
+    AddTriangleLine(TPS[0], TPS[2], 1);
+    AddTriangleLine(TPS[1], TPS[2], 2);
+
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    AddTriangleToLines();
+
+    BTS->GetNormalVector(BTS->NormalVector);
+
+    if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  || (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) || (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) ) {
+      BTS->NormalVector.Scale(-1);
+    };
+    cout << Verbose(1) << "New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
+    cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
+  } else { // else, yell and do nothing
+    cout << Verbose(1) << "This triangle consisting of ";
+    cout << *Opt_Candidate << ", ";
+    cout << *Line.endpoints[0]->node << " and ";
+    cout << *Line.endpoints[1]->node << " ";
+    cout << "is invalid!" << endl;
+    return false;
+  }
+
+  if ((TrianglesOnBoundaryCount % 10) == 0) {
+    sprintf(NumberName, "-%d", TriangleFilesWritten);
+    if (DoTecplotOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(TecplotSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+
+    if (DoRaster3DOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(Raster3DSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_raster3d_file(out, tempstream, this, mol);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoTecplotOutput || DoRaster3DOutput)
+      TriangleFilesWritten++;
+  }
+
+  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+  return true;
 };
 
@@ -2159,36 +2159,36 @@
  */
 bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c) {
-	LineMap::iterator TryAndError;
-	PointTestPair InsertPoint;
-	bool Present[3];
-
-	*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
-	TPS[0] = new class BoundaryPointSet(a);
-	TPS[1] = new class BoundaryPointSet(b);
-	TPS[2] = new class BoundaryPointSet(c);
-	for (int i=0;i<3;i++) { // check through all endpoints
-		InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i]));
-		Present[i] = !InsertPoint.second;
-		if (Present[i]) { // if new point was not present before, increase counter
-			delete TPS[i];
-			*out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl;
-			TPS[i] = (InsertPoint.first)->second;
-		}
-	}
-
-	// check lines
-	for (int i=0;i<3;i++)
-		if (Present[i])
-			for (int j=i;j<3;j++)
-				if (Present[j]) {
-					TryAndError = TPS[i]->lines.find(TPS[j]->node->nr);
-					if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) {
-						*out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl;
-						*out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
-						return false;
-					}
-				}
-	*out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
-	return true;
+  LineMap::iterator TryAndError;
+  PointTestPair InsertPoint;
+  bool Present[3];
+
+  *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+  TPS[0] = new class BoundaryPointSet(a);
+  TPS[1] = new class BoundaryPointSet(b);
+  TPS[2] = new class BoundaryPointSet(c);
+  for (int i=0;i<3;i++) { // check through all endpoints
+    InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i]));
+    Present[i] = !InsertPoint.second;
+    if (Present[i]) { // if new point was not present before, increase counter
+      delete TPS[i];
+      *out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl;
+      TPS[i] = (InsertPoint.first)->second;
+    }
+  }
+
+  // check lines
+  for (int i=0;i<3;i++)
+    if (Present[i])
+      for (int j=i;j<3;j++)
+        if (Present[j]) {
+          TryAndError = TPS[i]->lines.find(TPS[j]->node->nr);
+          if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) {
+            *out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl;
+            *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+            return false;
+          }
+        }
+  *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+  return true;
 };
 
@@ -2198,192 +2198,192 @@
     molecule* mol, double RADIUS)
 {
-	cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
-	int i;
-	Vector AngleCheck;
-	atom* Walker;
-	double norm = -1., angle;
-
-	// check if we only have one unique point yet ...
-	if (a != Candidate) {
-		cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
-		AngleCheck.CopyVector(&(Candidate->x));
-		AngleCheck.SubtractVector(&(a->x));
-		norm = AngleCheck.Norm();
-		// second point shall have smallest angle with respect to Oben vector
-		if (norm < RADIUS) {
-			angle = AngleCheck.Angle(&Oben);
-			if (angle < Storage[0]) {
-				//cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-				cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
-				Opt_Candidate = Candidate;
-				Storage[0] = AngleCheck.Angle(&Oben);
-				//cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-			} else {
-				cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
-			}
-		} else {
-			cout << "Refused due to Radius " << norm << endl;
-		}
-	}
-
-	// if not recursed to deeply, look at all its bonds
-	if (RecursionLevel < 7) {
-		for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) {
-			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-			if (Walker == Parent) // don't go back along the bond we came from
-				continue;
-			else
-				Find_second_point_for_Tesselation(a, Walker, Candidate, RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
-		}
-	}
-	cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
+  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
+  int i;
+  Vector AngleCheck;
+  atom* Walker;
+  double norm = -1., angle;
+
+  // check if we only have one unique point yet ...
+  if (a != Candidate) {
+    cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
+    AngleCheck.CopyVector(&(Candidate->x));
+    AngleCheck.SubtractVector(&(a->x));
+    norm = AngleCheck.Norm();
+    // second point shall have smallest angle with respect to Oben vector
+    if (norm < RADIUS) {
+      angle = AngleCheck.Angle(&Oben);
+      if (angle < Storage[0]) {
+        //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+        cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
+        Opt_Candidate = Candidate;
+        Storage[0] = AngleCheck.Angle(&Oben);
+        //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+      } else {
+        cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+      }
+    } else {
+      cout << "Refused due to Radius " << norm << endl;
+    }
+  }
+
+  // if not recursed to deeply, look at all its bonds
+  if (RecursionLevel < 7) {
+    for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) {
+      Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+      if (Walker == Parent) // don't go back along the bond we came from
+        continue;
+      else
+        Find_second_point_for_Tesselation(a, Walker, Candidate, RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+    }
+  }
+  cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
 };
 
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-	cout << Verbose(1) << "Begin of Find_starting_triangle\n";
-	int i = 0;
-	atom* Walker;
-	atom* FirstPoint;
-	atom* SecondPoint;
-	atom* max_index[NDIM];
-	double max_coordinate[NDIM];
-	Vector Oben;
-	Vector helper;
-	Vector Chord;
-	Vector CenterOfFirstLine;
-
-	Oben.Zero();
-
-	for (i = 0; i < 3; i++) {
-		max_index[i] = NULL;
-		max_coordinate[i] = -1;
-	}
-	cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
-
-	// 1. searching topmost atom with respect to each axis
-	Walker = mol->start;
-	while (Walker->next != mol->end) {
-		Walker = Walker->next;
-		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;
-			}
-		}
-	}
-
-	cout << Verbose(2) << "Found maximum coordinates: ";
-	for (int i=0;i<NDIM;i++)
-		cout << i << ": " << *max_index[i] << "\t";
-	cout << endl;
+  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+  int i = 0;
+  atom* Walker;
+  atom* FirstPoint;
+  atom* SecondPoint;
+  atom* max_index[NDIM];
+  double max_coordinate[NDIM];
+  Vector Oben;
+  Vector helper;
+  Vector Chord;
+  Vector CenterOfFirstLine;
+
+  Oben.Zero();
+
+  for (i = 0; i < 3; i++) {
+    max_index[i] = NULL;
+    max_coordinate[i] = -1;
+  }
+  cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
+
+  // 1. searching topmost atom with respect to each axis
+  Walker = mol->start;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    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;
+      }
+    }
+  }
+
+  cout << Verbose(2) << "Found maximum coordinates: ";
+  for (int i=0;i<NDIM;i++)
+    cout << i << ": " << *max_index[i] << "\t";
+  cout << endl;
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-	const int k = 1;
-	Oben.x[k] = 1.;
-	FirstPoint = max_index[k];
-
-	cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
-	double Storage[3];
-	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.
-	Storage[2] = 999999.;
-
-	Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
-	SecondPoint = Opt_Candidate;
-	cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-
-	helper.CopyVector(&(FirstPoint->x));
-	helper.SubtractVector(&(SecondPoint->x));
-	helper.Normalize();
-	Oben.ProjectOntoPlane(&helper);
-	Oben.Normalize();
-	helper.VectorProduct(&Oben);
-	Storage[0] = -2.; // This will indicate the quadrant.
-	Storage[1] = 9999999.; // This will be an angle looking for the third point.
-	Storage[2] = 9999999.;
-
-	Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-	Chord.SubtractVector(&(SecondPoint->x));
-	// Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-
-	cout << Verbose(2) << "Looking for third point candidates \n";
-	// look in one direction of baseline for initial candidate
-	Opt_Candidate = NULL;
-	CenterOfFirstLine.CopyVector(&Chord);
-	CenterOfFirstLine.Scale(0.5);
-	CenterOfFirstLine.AddVector(&(SecondPoint->x));
-
-	cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
-	Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
-	// look in other direction of baseline for possible better candidate
-	cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
-	Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
-	cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
-
-	// FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
-
-	// Finally, we only have to add the found points
-	AddTrianglePoint(FirstPoint, 0);
-	AddTrianglePoint(SecondPoint, 1);
-	AddTrianglePoint(Opt_Candidate, 2);
-	// ... and respective lines
-	AddTriangleLine(TPS[0], TPS[1], 0);
-	AddTriangleLine(TPS[1], TPS[2], 1);
-	AddTriangleLine(TPS[0], TPS[2], 2);
-	// ... and triangles to the Maps of the Tesselation class
-	BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-	AddTriangleToLines();
-	// ... and calculate its normal vector (with correct orientation)
-	Oben.Scale(-1.);
-	BTS->GetNormalVector(Oben);
-	cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
-	cout << Verbose(1) << "End of Find_starting_triangle\n";
+  const int k = 1;
+  Oben.x[k] = 1.;
+  FirstPoint = max_index[k];
+
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
+  double Storage[3];
+  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.
+  Storage[2] = 999999.;
+
+  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
+  SecondPoint = Opt_Candidate;
+  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+
+  helper.CopyVector(&(FirstPoint->x));
+  helper.SubtractVector(&(SecondPoint->x));
+  helper.Normalize();
+  Oben.ProjectOntoPlane(&helper);
+  Oben.Normalize();
+  helper.VectorProduct(&Oben);
+  Storage[0] = -2.; // This will indicate the quadrant.
+  Storage[1] = 9999999.; // This will be an angle looking for the third point.
+  Storage[2] = 9999999.;
+
+  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+  Chord.SubtractVector(&(SecondPoint->x));
+  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+
+  cout << Verbose(2) << "Looking for third point candidates \n";
+  // look in one direction of baseline for initial candidate
+  Opt_Candidate = NULL;
+  CenterOfFirstLine.CopyVector(&Chord);
+  CenterOfFirstLine.Scale(0.5);
+  CenterOfFirstLine.AddVector(&(SecondPoint->x));
+
+  cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
+  // look in other direction of baseline for possible better candidate
+  cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
+  cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+
+  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+
+  // Finally, we only have to add the found points
+  AddTrianglePoint(FirstPoint, 0);
+  AddTrianglePoint(SecondPoint, 1);
+  AddTrianglePoint(Opt_Candidate, 2);
+  // ... and respective lines
+  AddTriangleLine(TPS[0], TPS[1], 0);
+  AddTriangleLine(TPS[1], TPS[2], 1);
+  AddTriangleLine(TPS[0], TPS[2], 2);
+  // ... and triangles to the Maps of the Tesselation class
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTriangleToLines();
+  // ... and calculate its normal vector (with correct orientation)
+  Oben.Scale(-1.);
+  BTS->GetNormalVector(Oben);
+  cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
+  cout << Verbose(1) << "End of Find_starting_triangle\n";
 };
 
 void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol, const char *filename, const double RADIUS)
 {
-	int N = 0;
-	struct Tesselation *Tess = new Tesselation;
-	cout << Verbose(1) << "Entering search for non convex hull. " << endl;
-	cout << flush;
-	LineMap::iterator baseline;
-	cout << Verbose(0) << "Begin of Find_non_convex_border\n";
-	bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
-	bool failflag = false;
-	if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
-		mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
-
-	Tess->Find_starting_triangle(mol, RADIUS);
-
-	baseline = Tess->LinesOnBoundary.begin();
-	while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
-		if (baseline->second->TrianglesCount == 1) {
-			failflag = Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
-			flag = flag || failflag;
-			if (!failflag)
-				cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
-		} else {
-			cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
-		}
-		N++;
-		baseline++;
-		if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
-			baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-			flag = false;
-		}
-	}
-	if (failflag) {
-		cout << Verbose(1) << "Writing final tecplot file\n";
-		if (DoTecplotOutput)
-			write_tecplot_file(out, tecplot, Tess, mol, -1);
-		if (DoRaster3DOutput)
-			write_raster3d_file(out, tecplot, Tess, mol);
-	} else {
-		cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
-	}
-
-	cout << Verbose(0) << "End of Find_non_convex_border\n";
-	delete(Tess);
+  int N = 0;
+  struct Tesselation *Tess = new Tesselation;
+  cout << Verbose(1) << "Entering search for non convex hull. " << endl;
+  cout << flush;
+  LineMap::iterator baseline;
+  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
+  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
+  bool failflag = false;
+  if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
+    mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
+
+  Tess->Find_starting_triangle(mol, RADIUS);
+
+  baseline = Tess->LinesOnBoundary.begin();
+  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
+    if (baseline->second->TrianglesCount == 1) {
+      failflag = Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
+      flag = flag || failflag;
+      if (!failflag)
+        cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
+    } else {
+      cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
+    }
+    N++;
+    baseline++;
+    if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
+      baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      flag = false;
+    }
+  }
+  if (failflag) {
+    cout << Verbose(1) << "Writing final tecplot file\n";
+    if (DoTecplotOutput)
+      write_tecplot_file(out, tecplot, Tess, mol, -1);
+    if (DoRaster3DOutput)
+      write_raster3d_file(out, tecplot, Tess, mol);
+  } else {
+    cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
+  }
+
+  cout << Verbose(0) << "End of Find_non_convex_border\n";
+  delete(Tess);
 };
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/builder.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -562,16 +562,16 @@
       break;
     case 'd':
-    	cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
-    	cout << Verbose(0) << "Shall we rotate? [0/1]: ";
-    	cin >> Z;
-    	if ((Z >=0) && (Z <=1))
-    	  mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
-    	else
-    	  mol->PrincipalAxisSystem((ofstream *)&cout, false);
-    	break;
+      cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+      cout << Verbose(0) << "Shall we rotate? [0/1]: ";
+      cin >> Z;
+      if ((Z >=0) && (Z <=1))
+        mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
+      else
+        mol->PrincipalAxisSystem((ofstream *)&cout, false);
+      break;
     case 'e':
-    	cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-    	VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
-    	break;
+      cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+      VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
+      break;
     case 'f':
       mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
@@ -1020,33 +1020,33 @@
               }
               break;
-						case 'A':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (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 'N':
-							ExitFlag = 1;
-							if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
-							} else {
-								cout << Verbose(0) << "Evaluating npn-convex envelope.";
-								string TempName(argv[argptr+1]);
-								TempName.append(".r3d");
-								ofstream *output = new ofstream(TempName.c_str(), ios::trunc);
-								cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
-								Find_non_convex_border((ofstream *)&cout, output, mol, argv[argptr+1], atof(argv[argptr]));
-								output->close();
-								delete(output);
-								argptr+=2;
-							}
-							break;
+            case 'A':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (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 'N':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
+              } else {
+                cout << Verbose(0) << "Evaluating npn-convex envelope.";
+                string TempName(argv[argptr+1]);
+                TempName.append(".r3d");
+                ofstream *output = new ofstream(TempName.c_str(), ios::trunc);
+                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
+                Find_non_convex_border((ofstream *)&cout, output, mol, argv[argptr+1], atof(argv[argptr]));
+                output->close();
+                delete(output);
+                argptr+=2;
+              }
+              break;
             case 'T':
               ExitFlag = 1;
@@ -1369,5 +1369,5 @@
   char choice;  // menu choice char
   Vector x,y,z,n;  // coordinates for absolute point in cell volume
-	double *factor;	// unit factor if desired
+  double *factor;  // unit factor if desired
   bool valid; // flag if input was valid or not
   ifstream test;
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/config.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -1082,6 +1082,6 @@
 {
   bool result = true;
-	// bring MaxTypes up to date
-	mol->CountElements();
+  // bring MaxTypes up to date
+  mol->CountElements();
   ofstream *output = NULL;
   output = new ofstream(filename, ios::out);
Index: src/datacreator.cpp
===================================================================
--- src/datacreator.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/datacreator.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -335,6 +335,6 @@
     for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
       for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
-      	tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
-      	norm += tmp*tmp;
+        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
+        norm += tmp*tmp;
       }
     }
Index: src/helpers.cpp
===================================================================
--- src/helpers.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/helpers.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -73,5 +73,5 @@
   void *buffer = NULL;
   if (OldPointer == NULL)
-  	//cout << Verbose(0) << "ReAlloc impossible - old is NULL: " << output << endl;
+    //cout << Verbose(0) << "ReAlloc impossible - old is NULL: " << output << endl;
     buffer = (void *)malloc(size); // malloc
   else
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/molecules.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -53,5 +53,5 @@
   BondCount = 0;
   NoNonBonds = 0;
-	NoNonHydrogen = 0;
+  NoNonHydrogen = 0;
   NoCyclicBonds = 0;
   ListOfBondsPerAtom = NULL;
@@ -117,7 +117,7 @@
 {
   if (pointer != NULL) {
-  	atom *walker = new atom();
-  	walker->type = pointer->type;	// copy element of atom
-  	walker->x.CopyVector(&pointer->x); // copy coordination
+    atom *walker = new atom();
+    walker->type = pointer->type;  // copy element of atom
+    walker->x.CopyVector(&pointer->x); // copy coordination
     walker->v.CopyVector(&pointer->v); // copy velocity
     walker->FixedIon = pointer->FixedIon;
@@ -186,5 +186,5 @@
   double *matrix;
 
-//	*out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
+//  *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
   InBondvector.CopyVector(&TopReplacement->x);
@@ -427,5 +427,5 @@
   }
 
-//	*out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
+//  *out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
   return AllWentWell;
 };
@@ -753,10 +753,10 @@
 Vector * molecule::DetermineCenterOfGravity(ofstream *out)
 {
-	atom *ptr = start->next;  // start at first in list
-	Vector *a = new Vector();
-	Vector tmp;
+  atom *ptr = start->next;  // start at first in list
+  Vector *a = new Vector();
+  Vector tmp;
   double Num = 0;
 
-	a->Zero();
+  a->Zero();
 
   if (ptr != end) {   //list not empty?
@@ -909,68 +909,68 @@
 void molecule::PrincipalAxisSystem(ofstream *out, bool DoRotate)
 {
-	atom *ptr = start;  // start at first in list
-	double InertiaTensor[NDIM*NDIM];
-	Vector *CenterOfGravity = DetermineCenterOfGravity(out);
-
-	CenterGravity(out, CenterOfGravity);
-
-	// reset inertia tensor
-	for(int i=0;i<NDIM*NDIM;i++)
-		InertiaTensor[i] = 0.;
-
-	// sum up inertia tensor
-	while (ptr->next != end) {
-		ptr = ptr->next;
-		Vector x;
-		x.CopyVector(&ptr->x);
-		//x.SubtractVector(CenterOfGravity);
-		InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-		InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-		InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-		InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-		InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-		InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-		InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-		InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-		InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
-	}
-	// print InertiaTensor for debugging
-	*out << "The inertia tensor is:" << endl;
-	for(int i=0;i<NDIM;i++) {
-	  for(int j=0;j<NDIM;j++)
-	    *out << InertiaTensor[i*NDIM+j] << " ";
-	  *out << endl;
-	}
-	*out << endl;
-
-	// diagonalize to determine principal axis system
-	gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
-	gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
-	gsl_vector *eval = gsl_vector_alloc(NDIM);
-	gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
-	gsl_eigen_symmv(&m.matrix, eval, evec, T);
-	gsl_eigen_symmv_free(T);
-	gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
-
-	for(int i=0;i<NDIM;i++) {
-		*out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
-		*out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
-	}
-
-	// check whether we rotate or not
-	if (DoRotate) {
-	  *out << Verbose(1) << "Transforming molecule into PAS ... ";
-	  // the eigenvectors specify the transformation matrix
-	  ptr = start;
-	  while (ptr->next != end) {
-	    ptr = ptr->next;
-	    for (int j=0;j<MDSteps;j++)
-	      Trajectories[ptr].R.at(j).MatrixMultiplication(evec->data);
-	    ptr->x.MatrixMultiplication(evec->data);
-	  }
-	  *out << "done." << endl;
-
-	  // summing anew for debugging (resulting matrix has to be diagonal!)
-	  // reset inertia tensor
+  atom *ptr = start;  // start at first in list
+  double InertiaTensor[NDIM*NDIM];
+  Vector *CenterOfGravity = DetermineCenterOfGravity(out);
+
+  CenterGravity(out, CenterOfGravity);
+
+  // reset inertia tensor
+  for(int i=0;i<NDIM*NDIM;i++)
+    InertiaTensor[i] = 0.;
+
+  // sum up inertia tensor
+  while (ptr->next != end) {
+    ptr = ptr->next;
+    Vector x;
+    x.CopyVector(&ptr->x);
+    //x.SubtractVector(CenterOfGravity);
+    InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
+    InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
+    InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
+    InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
+    InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
+    InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
+    InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
+    InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
+    InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+  }
+  // print InertiaTensor for debugging
+  *out << "The inertia tensor is:" << endl;
+  for(int i=0;i<NDIM;i++) {
+    for(int j=0;j<NDIM;j++)
+      *out << InertiaTensor[i*NDIM+j] << " ";
+    *out << endl;
+  }
+  *out << endl;
+
+  // diagonalize to determine principal axis system
+  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
+  gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
+  gsl_vector *eval = gsl_vector_alloc(NDIM);
+  gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
+  gsl_eigen_symmv(&m.matrix, eval, evec, T);
+  gsl_eigen_symmv_free(T);
+  gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
+
+  for(int i=0;i<NDIM;i++) {
+    *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
+    *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
+  }
+
+  // check whether we rotate or not
+  if (DoRotate) {
+    *out << Verbose(1) << "Transforming molecule into PAS ... ";
+    // the eigenvectors specify the transformation matrix
+    ptr = start;
+    while (ptr->next != end) {
+      ptr = ptr->next;
+      for (int j=0;j<MDSteps;j++)
+        Trajectories[ptr].R.at(j).MatrixMultiplication(evec->data);
+      ptr->x.MatrixMultiplication(evec->data);
+    }
+    *out << "done." << endl;
+
+    // summing anew for debugging (resulting matrix has to be diagonal!)
+    // reset inertia tensor
     for(int i=0;i<NDIM*NDIM;i++)
       InertiaTensor[i] = 0.;
@@ -1001,10 +1001,10 @@
     }
     *out << endl;
-	}
-
-	// free everything
-	delete(CenterOfGravity);
-	gsl_vector_free(eval);
-	gsl_matrix_free(evec);
+  }
+
+  // free everything
+  delete(CenterOfGravity);
+  gsl_vector_free(eval);
+  gsl_matrix_free(evec);
 };
 
@@ -2027,5 +2027,5 @@
     runner = elemente->start;
     while (runner->next != elemente->end) { // go through every element
-  		runner = runner->next;
+      runner = runner->next;
       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
         ElementNo++;
@@ -2121,5 +2121,5 @@
 bool molecule::Checkout(ofstream *out)  const
 {
-	return elemente->Checkout(out, ElementsInMolecule);
+  return elemente->Checkout(out, ElementsInMolecule);
 };
 
@@ -2174,5 +2174,5 @@
     runner = elemente->start;
     while (runner->next != elemente->end) { // go through every element
-  		runner = runner->next;
+      runner = runner->next;
       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
         ElementNo++;
@@ -2196,6 +2196,6 @@
 void molecule::CountAtoms(ofstream *out)
 {
-	int i = 0;
-	atom *Walker = start;
+  int i = 0;
+  atom *Walker = start;
   while (Walker->next != end) {
     Walker = Walker->next;
@@ -2204,24 +2204,24 @@
   if ((AtomCount == 0) || (i != AtomCount)) {
     *out << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
-  	AtomCount = i;
+    AtomCount = i;
 
     // count NonHydrogen atoms and give each atom a unique name
     if (AtomCount != 0) {
-	    i=0;
-	    NoNonHydrogen = 0;
+      i=0;
+      NoNonHydrogen = 0;
       Walker = start;
-	    while (Walker->next != end) {
-	      Walker = Walker->next;
-	      Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-	      if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
-	        NoNonHydrogen++;
-	      Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
-	      Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
-	      sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
+      while (Walker->next != end) {
+        Walker = Walker->next;
+        Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
+        if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
+          NoNonHydrogen++;
+        Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
+        Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
+        sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
         *out << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
-	      i++;
-	    }
+        i++;
+      }
     } else
-    	*out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
+      *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
   }
 };
@@ -2231,8 +2231,8 @@
 void molecule::CountElements()
 {
-	int i = 0;
+  int i = 0;
   for(i=MAX_ELEMENTS;i--;)
-  	ElementsInMolecule[i] = 0;
-	ElementCount = 0;
+    ElementsInMolecule[i] = 0;
+  ElementCount = 0;
 
   atom *walker = start;
@@ -2243,5 +2243,5 @@
   }
   for(i=MAX_ELEMENTS;i--;)
-  	ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
+    ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
 };
 
@@ -2334,42 +2334,42 @@
 {
 
-	// 1 We will parse bonds out of the dbond file created by tremolo.
-			int atom1, atom2, temp;
-			atom *Walker, *OtherWalker;
-
-	        if (!input)
-	        {
-	        	cout << Verbose(1) << "Opening silica failed \n";
-	        };
-
-			*input >> ws >> atom1;
-			*input >> ws >> atom2;
-	        cout << Verbose(1) << "Scanning file\n";
-	        while (!input->eof()) // Check whether we read everything already
-	        {
-				*input >> ws >> atom1;
-				*input >> ws >> atom2;
-	        	if(atom2<atom1) //Sort indices of atoms in order
-	        	{
-	        		temp=atom1;
-	        		atom1=atom2;
-	        		atom2=temp;
-	        	};
-
-	        	Walker=start;
-	        	while(Walker-> nr != atom1) // Find atom corresponding to first index
-	        	{
-	        		Walker = Walker->next;
-	        	};
-	        	OtherWalker = Walker->next;
-	        	while(OtherWalker->nr != atom2) // Find atom corresponding to second index
-	        	{
-	        		OtherWalker= OtherWalker->next;
-	        	};
-	        	AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
-	        	
-	        }
-
-	        CreateListOfBondsPerAtom(out);
+  // 1 We will parse bonds out of the dbond file created by tremolo.
+      int atom1, atom2, temp;
+      atom *Walker, *OtherWalker;
+
+          if (!input)
+          {
+            cout << Verbose(1) << "Opening silica failed \n";
+          };
+
+      *input >> ws >> atom1;
+      *input >> ws >> atom2;
+          cout << Verbose(1) << "Scanning file\n";
+          while (!input->eof()) // Check whether we read everything already
+          {
+        *input >> ws >> atom1;
+        *input >> ws >> atom2;
+            if(atom2<atom1) //Sort indices of atoms in order
+            {
+              temp=atom1;
+              atom1=atom2;
+              atom2=temp;
+            };
+
+            Walker=start;
+            while(Walker-> nr != atom1) // Find atom corresponding to first index
+            {
+              Walker = Walker->next;
+            };
+            OtherWalker = Walker->next;
+            while(OtherWalker->nr != atom2) // Find atom corresponding to second index
+            {
+              OtherWalker= OtherWalker->next;
+            };
+            AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
+            
+          }
+
+          CreateListOfBondsPerAtom(out);
 
 };
@@ -2520,9 +2520,9 @@
     // 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) {
+    if (BondCount != 0) {
       NoCyclicBonds = 0;
-		  *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
-		  do {
-		    No = 0; // No acts as breakup flag (if 1 we still continue)
+      *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
+      do {
+        No = 0; // No acts as breakup flag (if 1 we still continue)
         Walker = start;
         while (Walker->next != end) { // go through every atom
@@ -2538,17 +2538,17 @@
             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
+              // 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 " << *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
+              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;
-	              }
-	            }
-	          }
+                  *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
+                }
+              }
+            }
             if ((Candidate != NULL) && (CandidateBondNo != -1)) {
               ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
@@ -2558,21 +2558,21 @@
               FalseBondDegree++;
           }
-		    }
-		  } while (No);
-		*out << " done." << endl;
-		} else
-			*out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
-	  *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
-		
-	  // output bonds for debugging (if bond chain list was correctly installed)
-	  *out << Verbose(1) << endl << "From contents of bond chain list:";
-	  bond *Binder = first;
+        }
+      } while (No);
+    *out << " done." << endl;
+    } else
+      *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
+    *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
+    
+    // output bonds for debugging (if bond chain list was correctly installed)
+    *out << Verbose(1) << endl << "From contents of bond chain list:";
+    bond *Binder = first;
     while(Binder->next != last) {
       Binder = Binder->next;
-			*out << *Binder << "\t" << endl;
+      *out << *Binder << "\t" << endl;
     }
     *out << endl;
   } else
-  	*out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
+    *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
   Free((void **)&matrix, "molecule::CreateAdjacencyList: *matrix");
@@ -3600,5 +3600,5 @@
 int molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
 {
-	MoleculeListClass *BondFragments = NULL;
+  MoleculeListClass *BondFragments = NULL;
   int *SortIndex = NULL;
   int *MinimumRingSize = new int[AtomCount];
@@ -3806,12 +3806,12 @@
     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]);
+      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;
-	        break;
-  	    }
-    	}
+          break;
+        }
+      }
     Binder = ReferenceStack->PopFirst();  // loop the stack for next item
     *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
@@ -4308,5 +4308,5 @@
   for (int i=0;i<AtomCount;i++) { // reset all atom labels
     // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
-		Labels[i] = -1;
+    Labels[i] = -1;
     SonList[i] = NULL;
     PredecessorList[i] = NULL;
@@ -4316,5 +4316,5 @@
   for (int i=0;i<BondCount;i++)
     ColorEdgeList[i] = white;
-	RootStack->ClearStack();	// clearstack and push first atom if exists
+  RootStack->ClearStack();  // clearstack and push first atom if exists
   TouchedStack->ClearStack();
   Walker = start->next;
@@ -4335,6 +4335,6 @@
   ///// OUTER LOOP ////////////
   while (!RootStack->IsEmpty()) {
-  	// get new root vertex from atom stack
-  	Root = RootStack->PopFirst();
+    // get new root vertex from atom stack
+    Root = RootStack->PopFirst();
     ShortestPathList[Root->nr] = 0;
     if (Labels[Root->nr] == -1)
@@ -4344,6 +4344,6 @@
     *out << Verbose(0) << "Root for this loop is: " << Root->Name << ".\n";
 
-		// clear snake stack
-	  SnakeStack->ClearStack();
+    // clear snake stack
+    SnakeStack->ClearStack();
     //SnakeStack->TestImplementation(out, start->next);
 
@@ -4352,13 +4352,13 @@
     // - what about cyclic bonds?
     Walker = Root;
-  	do {
+    do {
       *out << Verbose(1) << "Current Walker is: " << Walker->Name;
       // initial setting of the new Walker: label, color, shortest path and put on stacks
-  		if (Labels[Walker->nr] == -1)	{	// give atom a unique, monotonely increasing number
-  			Labels[Walker->nr] = RunningIndex++;
+      if (Labels[Walker->nr] == -1)  {  // give atom a unique, monotonely increasing number
+        Labels[Walker->nr] = RunningIndex++;
         RootStack->Push(Walker);
       }
       *out << ", has label " << Labels[Walker->nr];
-  		if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) {	// color it if newly discovered and push on stacks (and if within reach!)
+      if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) {  // color it if newly discovered and push on stacks (and if within reach!)
         if ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white)) {
           // Binder ought to be set still from last neighbour search
@@ -4374,9 +4374,9 @@
           ColorVertexList[Walker->nr] = darkgray; // mark as dark gray of on snake stack
         }
-  		}
+      }
       *out << ", SP of " << ShortestPathList[Walker->nr]  << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
 
       // then check the stack for a newly stumbled upon fragment
-  		if (SnakeStack->ItemCount() == Order) { // is stack full?
+      if (SnakeStack->ItemCount() == Order) { // is stack full?
         // store the fragment if it is one and get a removal candidate
         Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
@@ -4391,5 +4391,5 @@
           }
         }
-  		} else
+      } else
         Removal = NULL;
 
@@ -4398,15 +4398,15 @@
       if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {  // don't look, if a new walker has been set above
         *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
-  			OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
+        OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
         if (ShortestPathList[Walker->nr] < Order) {
-    			for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
+          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
             Binder = ListOfBondsPerAtom[Walker->nr][i];
             *out << Verbose(2) << "Current bond is " << *Binder << ": ";
-    				OtherAtom = Binder->GetOtherAtom(Walker);
+            OtherAtom = Binder->GetOtherAtom(Walker);
             if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
               *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
               //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
             } else { // otherwise check its colour and element
-      				if (
+              if (
 #ifdef ADDHYDROGEN
               (OtherAtom->type->Z != 1) &&
@@ -4430,5 +4430,5 @@
               }
             }
-    			}
+          }
         } else {  // means we have stepped beyond the horizon: Return!
           Walker = PredecessorList[Walker->nr];
@@ -4436,11 +4436,11 @@
           *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
         }
-  			if (Walker != OtherAtom) {	// if no white neighbours anymore, color it black
+        if (Walker != OtherAtom) {  // if no white neighbours anymore, color it black
           *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
-  				ColorVertexList[Walker->nr] = black;
-  				Walker = PredecessorList[Walker->nr];
-  			}
-      }
-  	} while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
+          ColorVertexList[Walker->nr] = black;
+          Walker = PredecessorList[Walker->nr];
+        }
+      }
+    } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
     *out << Verbose(2) << "Inner Looping is finished." << endl;
 
@@ -4562,5 +4562,5 @@
         bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
         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
+          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][i] << ", checking on " << *OtherWalker << "." << endl;
           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
@@ -4584,9 +4584,9 @@
         if (SubOrder > 1) {    // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
           // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
-          SP = RootDistance+1;	// this is the next level
+          SP = RootDistance+1;  // this is the next level
           // first count the members in the subset
           SubSetDimension = 0;
-          Binder = FragmentSearch->BondsPerSPList[2*SP];		// start node for this level
-          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {		// compare to end node of this level
+          Binder = FragmentSearch->BondsPerSPList[2*SP];    // start node for this level
+          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {    // compare to end node of this level
             Binder = Binder->next;
             for (int k=TouchedIndex;k--;) {
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/parser.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -786,5 +786,5 @@
           }
           if (k != -1) {
-          	//cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
+            //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
             Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
           }
@@ -855,8 +855,8 @@
                   }
                   if (Order == SubOrder) { // equal order is always copy from Energies
-                  	//cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
                   } else {
-                  	//cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
                   }
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/vector.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -327,5 +327,5 @@
   }
   ost << ")";
-	return ost;
+  return ost;
 };
 
@@ -517,5 +517,5 @@
 bool Vector::MakeNormalVector(const Vector *y1)
 {
-	bool result = false;
+  bool result = false;
   Vector x1;
   x1.CopyVector(y1);
@@ -523,5 +523,5 @@
   SubtractVector(&x1);
   for (int i=NDIM;i--;)
-	  result = result || (fabs(x[i]) > MYEPSILON);
+    result = result || (fabs(x[i]) > MYEPSILON);
 
   return result;
@@ -597,14 +597,14 @@
 bool Vector::LSQdistance(Vector **vectors, int num)
 {
-	int j;
-
-	for (j=0;j<num;j++) {
-		cout << Verbose(1) << j << "th atom's vector: ";
-		(vectors[j])->Output((ofstream *)&cout);
-		cout << endl;
-	}
+  int j;
+
+  for (j=0;j<num;j++) {
+    cout << Verbose(1) << j << "th atom's vector: ";
+    (vectors[j])->Output((ofstream *)&cout);
+    cout << endl;
+  }
 
   int np = 3;
-	struct LSQ_params par;
+  struct LSQ_params par;
 
    const gsl_multimin_fminimizer_type *T =
@@ -627,8 +627,8 @@
    /* Starting point */
    par.vectors = vectors;
-	 par.num = num;
-
-	 for (i=NDIM;i--;)
-		gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
+   par.num = num;
+
+   for (i=NDIM;i--;)
+    gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
 
    /* Initialize method and iterate */
@@ -671,5 +671,5 @@
    gsl_multimin_fminimizer_free (s);
 
-	return true;
+  return true;
 };
 
Index: src/verbose.cpp
===================================================================
--- src/verbose.cpp	(revision 986c809f30e6539ae660442a47de63b08aa35574)
+++ src/verbose.cpp	(revision ce5ac38e73b73390de3143a6985456a9e35495f9)
@@ -22,5 +22,5 @@
 ostream& operator<<(ostream& ost,const Verbose& m)
 {
-	return m.print(ost);
+  return m.print(ost);
 };
 
@@ -33,12 +33,12 @@
 ostream& Binary::print (ostream &ost) const
 {
-	int bits = 1, counter = 1;
-	while ((bits = 1 << counter) < BinaryNumber)
+  int bits = 1, counter = 1;
+  while ((bits = 1 << counter) < BinaryNumber)
     counter++;
   for (int i=0;i<counter-1;i++) {
-  	if ((BinaryNumber & (1 << i)) == 0)
-  		ost.put('0');
- 		else
-  		ost.put('1');
+    if ((BinaryNumber & (1 << i)) == 0)
+      ost.put('0');
+     else
+      ost.put('1');
   }
   ost.put('b');
@@ -56,4 +56,4 @@
 ostream& operator<<(ostream& ost,const Binary& m)
 {
-	return m.print(ost);
+  return m.print(ost);
 };
