Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 233e33e9480b84b60d83d44696ee46721204ee12)
+++ src/molecules.cpp	(revision 49de64f43390d73584f072908bf52bbacde73c26)
@@ -1381,6 +1381,14 @@
 	
 	// 1. calculate center of gravity
+	*out << endl;
 	vector *CenterOfGravity = DetermineCenterOfGravity(out);
 	
+	// 2. translate all points into CoG
+	*out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+	Walker=start;
+	while (Walker->next != end) {
+	  Walker = Walker->next;
+	  Walker->x.Translate(CenterOfGravity);
+	}
 	// 2. calculate distance of each atom and sort into hash table
 //	map<double, int> DistanceSet;
@@ -1393,11 +1401,25 @@
 	
 	// 3. Find all points on the boundary
+  *out << Verbose(1) << "Finding all boundary points." << endl;
 	Boundaries BoundaryPoints[NDIM];	// first is alpha, second is (r, nr)
 	BoundariesTestPair BoundaryTestPair;
-	vector AxisVector;
+	vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+	double radius, angle;
 	// 3a. Go through every axis
 	for (int axis=0; axis<NDIM; axis++)  {
 		AxisVector.Zero();
+    AngleReferenceVector.Zero();
+    AngleReferenceNormalVector.Zero();
 		AxisVector.x[axis] = 1.;
+		AngleReferenceVector.x[(axis+1)%NDIM] = 1.;
+		AngleReferenceNormalVector.x[(axis+2)%NDIM] = 1.;
+//		*out << Verbose(1) << "Axisvector is ";
+//		AxisVector.Output(out);
+//		*out << " and AngleReferenceVector is ";
+//		AngleReferenceVector.Output(out);
+//		*out << "." << endl;
+//    *out << " and AngleReferenceNormalVector is ";
+//    AngleReferenceNormalVector.Output(out);
+//    *out << "." << endl;
 		// 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
 		Walker = start;
@@ -1406,9 +1428,24 @@
 			vector ProjectedVector;
 			ProjectedVector.CopyVector(&Walker->x);
-			ProjectedVector.Scale(-1.*Walker->x.Projection(&AxisVector));
-			ProjectedVector.AddVector(&Walker->x);	// subtract projection from vector
-			BoundaryTestPair = BoundaryPoints[axis].insert( BoundariesPair (ProjectedVector.Angle(&AxisVector), DistanceNrPair (ProjectedVector.Norm(), Walker) ) );
+			ProjectedVector.ProjectOntoPlane(&AxisVector);
+		  // correct for negative side
+		  //if (Projection(y) < 0)
+		    //angle = 2.*M_PI - angle;
+			radius = ProjectedVector.Norm();
+			if (fabs(radius) > MYEPSILON)
+			  angle = ProjectedVector.Angle(&AngleReferenceVector);
+			else
+			  angle = 0.;  // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+			  
+			//*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+			if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+			  angle = 2.*M_PI - angle;
+			}
+			*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
+			ProjectedVector.Output(out);
+			*out << endl;
+			BoundaryTestPair = BoundaryPoints[axis].insert( BoundariesPair (angle, DistanceNrPair (radius, Walker) ) );
 			if (BoundaryTestPair.second) { // successfully inserted
-			} else { // same point exists, check first r, then distance of original vectors to centers of gravity
+			} else { // same point exists, check first r, then distance of original vectors to center of gravity
 				*out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
 				*out << Verbose(2) << "Present vector: ";
@@ -1424,5 +1461,5 @@
 					*out << Verbose(2) << "Keeping new vector." << endl;
 				} else if (tmp == BoundaryTestPair.first->second.first) {
-					if (BoundaryTestPair.first->second.second->x.Distance(CenterOfGravity) < Walker->x.Distance(CenterOfGravity)) {
+					if (BoundaryTestPair.first->second.second->x.ScalarProduct(&BoundaryTestPair.first->second.second->x) < Walker->x.ScalarProduct(&Walker->x)) { // Norm() does a sqrt, which makes it a lot slower
 						BoundaryTestPair.first->second.second = Walker;
 						*out << Verbose(2) << "Keeping new vector." << endl;
@@ -1435,7 +1472,21 @@
 			}
 		}
-		// 3c. throw out points whose distance is less than the mean of left and right neighbours
+		// printing all inserted for debugging
+		{
+      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      int i=0;
+      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        if (runner != BoundaryPoints[axis].begin())
+          *out << ", " << i << ": " << *runner->second.second;
+        else
+          *out << i << ": " << *runner->second.second;
+        i++;
+      }
+      *out << endl;
+		}
+    // 3c. throw out points whose distance is less than the mean of left and right neighbours
 		bool flag = false;
 		do { // do as long as we still throw one out per round
+		  *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
 			flag = false;
 			Boundaries::iterator left = BoundaryPoints[axis].end();
@@ -1455,8 +1506,52 @@
 				}
 				// check distance
-				if (runner->second.first < (left->second.first + right->second.first)/2. ) {
-					// throw out point
-					BoundaryPoints[axis].erase(runner);
-					flag = true;
+				
+				// construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+				{
+          vector SideA, SideB, SideC, SideH;
+          SideA.CopyVector(&left->second.second->x);
+          SideA.ProjectOntoPlane(&AxisVector);
+//          *out << "SideA: ";
+//          SideA.Output(out);
+//          *out << endl;
+          
+          SideB.CopyVector(&right->second.second->x);
+          SideB.ProjectOntoPlane(&AxisVector);
+//          *out << "SideB: ";
+//          SideB.Output(out);
+//          *out << endl;
+          
+          SideC.CopyVector(&left->second.second->x);
+          SideC.SubtractVector(&right->second.second->x);
+          SideC.ProjectOntoPlane(&AxisVector);
+//          *out << "SideC: ";
+//          SideC.Output(out);
+//          *out << endl;
+
+          SideH.CopyVector(&runner->second.second->x);
+          SideH.ProjectOntoPlane(&AxisVector);
+//          *out << "SideH: ";
+//          SideH.Output(out);
+//          *out << endl;
+          
+          // calculate each length
+          double a = SideA.Norm();
+          //double b = SideB.Norm();
+          //double c = SideC.Norm();
+          double h = SideH.Norm();
+          // calculate the angles
+          double alpha = SideA.Angle(&SideH);
+  				double beta = SideA.Angle(&SideC); 
+  				double gamma = SideB.Angle(&SideH);
+  				double delta = SideC.Angle(&SideH);
+  				double MinDistance = a * sin(beta)/(sin(delta)) * (((alpha < M_PI/2.) || (gamma < M_PI/2.)) ? 1. : -1.);
+//  				*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl; 
+  				*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+  				if ((fabs(h/fabs(h) - MinDistance/fabs(MinDistance)) < MYEPSILON) && (h <  MinDistance)) {
+  					// throw out point
+  				  *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+  					BoundaryPoints[axis].erase(runner);
+  					flag = true;
+  				}
 				}
 			}
@@ -1472,20 +1567,27 @@
 		}
 	}
-	int BoundaryPointCount = 0;
-	Walker = start;
-	while (Walker->next != end) {
-		Walker = Walker->next;
-		if (AtomList[Walker->nr] < NDIM) {	// enter all neighbouring points
-			BoundaryPointCount++;
-		}
-	} 
-	
+
 	// 3e. Points no more in the total list, have to be thrown out of each axis lists too!
 	for (int axis=0; axis<NDIM; axis++)  {
 		for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-			if (AtomList[runner->second.second->nr] < NDIM)
+			if (AtomList[runner->second.second->nr] < 1) {
+			  *out << Verbose(1) << "Throwing out " << *runner->second.second << " in axial projection of axis " << axis << "." << endl;
 				BoundaryPoints[axis].erase(runner);
+			}
 		}
 	}
+	
+	// listing boundary points
+	*out << Verbose(1) << "The following atoms are on the boundary:";
+  int BoundaryPointCount = 0;
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if (AtomList[Walker->nr] > 0) {
+      BoundaryPointCount ++;
+      *out << " " << Walker->Name;
+    }
+  }
+  *out << endl;
 
 	*out << Verbose(2) << "I found " << BoundaryPointCount << " points on the convex boundary." << endl;
@@ -1493,100 +1595,130 @@
 	
 	// 4. Create each possible line, number them and put them into a list (3 * AtomCount possible)
-	struct lines {
-		atom *x[2];
-		int nr;
-	} LinesList[3 * AtomCount];
-	
-	// initialise reference storage (needed lateron)
-	int **AtomLines = new int *[AtomCount];
-	Walker = start;
-	while (Walker->next != end) {	// go through all points
-		Walker = Walker->next;
-		if (AtomList[Walker->nr] == NDIM) {	// if this is a boundary point
-			AtomLines[Walker->nr] = new int[NDIM];
-			for(int axis=0;axis<NDIM;axis++) {
-				AtomLines[Walker->nr][2*axis+0] = -1;
-				AtomLines[Walker->nr][2*axis+1] = -1;
-			}
-		}
-	}
+	LinesMap AllTriangleLines;
 	
 	// then create each line
+  *out << Verbose(1) << "Creating lines between adjacent boundary points." << endl;
+  atom *Runner = NULL;
   int LineCount = 0;
-	for(int axis=0;axis<NDIM;axis++) {	// go through every axis
-	  for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-	    // get left and right neighbours
-	    Boundaries::iterator NN[2] = runner;
-	    if (NN[0] == BoundaryPoints[axis].begin())
-	      NN[0] = BoundaryPoints[axis].end();
-	    NN[0]--;
-	    NN[1] = runner;
-      NN[1]++;
-      if (NN[1] == BoundaryPoints[axis].end())
-        NN[1] = BoundaryPoints[axis].begin();
-      // check those neighbours for their atomic index
-      for (int k=0;k<2;k++) {
-  			if (runner->second.second->nr > NN[k]->second.second->nr) { // check left neighbour on this projected axis
-  				LinesList[LineCount].x[0] = runner->second.second;
-  				LinesList[LineCount].x[1] = NN[k]->second.second;
-  				if (AtomLines[LinesList[LineCount].x[0]->nr][axis] == -1) {
-  					AtomLines[LinesList[LineCount].x[0]->nr][axis] = LineCount;
-  				} else if (AtomLines[LinesList[LineCount].x[1]->nr][axis] == -1) {
-  					AtomLines[LinesList[LineCount].x[1]->nr][axis] = LineCount;
-  				} else { 
-  					*out << Verbose(2) << "Could not store line number " << LineCount << "." << endl;
-  				}
-  				LineCount++;
-  			}
-      }
-		}
-	}
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if (AtomList[Walker->nr] > 0) { // boundary point!
+      *out << Verbose(1) << "Current Walker is " << *Walker << "." << endl;
+      // make a list of all neighbours
+      DistanceMap Distances;
+      double distance = 0;
+      Runner = start;
+      while (Runner->next != end) {
+        Runner = Runner->next;
+        if ((AtomList[Runner->nr] > 0) && (Runner != Walker)) { // don't mess with ourselves!
+          distance = Walker->x.Distance(&Runner->x);    // note this is squared distance
+          Distances.insert ( DistanceNrPair( distance, Runner) );
+          *out << Verbose(2) << "Inserted distance " << distance << " to " << *Runner << " successfully." << endl;
+        }
+      } // end of distance filling while loop
+      // take 3NNs to draw a line in between
+      int Counter = 0;  // counts up to three inserted lines
+      LinesTestPair LinesTest;
+      for (DistanceMap::iterator spanner = Distances.begin(); spanner != Distances.end(); spanner++) {
+        *out << Verbose(1) << "Current distance to insert is " << spanner->first << " to " << spanner->second+1 << "." << endl;
+        LinesTest = AllTriangleLines.insert ( LinesPair( (Walker->nr*AtomCount)+spanner->second->nr, LineCount) );
+        if (!LinesTest.second) {
+          *out << Verbose(2) << "Insertion of line between " << Walker->nr+1 << " and " << spanner->second->nr+1 << " failed, already present line nr. " << LinesTest.first->second << "." << endl; 
+        } else {
+          LineCount++;
+          Counter++;
+          *out << Verbose(2) << "Insertion of line between " << (LinesTest.first->first/AtomCount)+1 << " and " << (LinesTest.first->first%AtomCount)+1 << " successful, inserting mirrored line and increased counter to " << Counter << "." << endl; 
+          LinesTest = AllTriangleLines.insert ( LinesPair( Walker->nr+(spanner->second->nr*AtomCount), LineCount) );
+          if (!LinesTest.second) {
+            *out << Verbose(2) << "Insertion of mirror line between " << Walker->nr << " and " << spanner->second->nr << " failed, already present line nr. " << LinesTest.first->second << "." << endl;
+            *out << Verbose(0) << "SERIOUS ERROR!!!" << endl;
+          }
+        }
+        if (Counter >= 3) { // stop after three lines (if not run out of possible NNs already)
+          *out << Verbose(2) << "Ok, three lines inserted from current walker " << *Walker << ", stopping." << endl;
+          break;
+        }
+      }
+    } // end of if boundary pointy
+  } // end of loop through all atoms
+
+	// listing created lines
+	*out << Verbose(2) << "The following lines were created:";
+	for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++)
+	  if ((liner->first/AtomCount) < (liner->first%AtomCount))
+	    *out << " " << (liner->first/AtomCount)+1 << "<->" << (liner->first%AtomCount)+1; 
+	*out << endl;
 	*out << Verbose(2) << "I created " << LineCount << " edges." << endl;
 	
 	// 5. Create every possible triangle from the lines (numbering of each line must be i < j < k)
 	struct triangles {
-		atom *x1;
-		atom *x2;
-		atom *x3;
+		atom *x[3];
 		int nr;
 	} TriangleList[2 * 3 * AtomCount];	// each line can be part in at most 2 triangles
+
+  *out << Verbose(1) << "Creating triangles out of these lines." << endl;
 	int TriangleCount = 0;
-	for (int i=0; i<LineCount; i++) {	// go through every line
+  LinetoAtomPair InsertionPair = LinetoAtomPair(-1, NULL);
+  int endpoints[3];
+  for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++) {	// go through every line
 		LinetoAtomMap LinetoAtom;
 		// we have two points, check the lines at either end, whether they share a common endpoint
-		for (int endpoint=0;endpoint<2;endpoint++) {
-			Walker = LinesList[LineCount].x[endpoint];		
-			for (int axis=0;axis<NDIM;axis++) {
-				LinetoAtomTestPair LinetoAtomTest;
-				LinetoAtomPair InsertionPair = LinetoAtomPair(-1, NULL);
-				struct lines *TempLine = &LinesList[ AtomLines[Walker->nr][axis] ];
-				if (AtomLines[Walker->nr][axis] != LineCount) {
-					if (TempLine->x[0] != Walker) {
-						InsertionPair = LinetoAtomPair(TempLine->x[0]->nr, TempLine->x[0]);
-					} else if (TempLine->x[1] != Walker) {
-						InsertionPair = LinetoAtomPair(TempLine->x[1]->nr, TempLine->x[0]);
-					} else {
-						*out << Verbose(2) << "Atom " << *Walker << "is both the end of line " << AtomLines[Walker->nr][axis] << "." << endl;
-					}
-					if (InsertionPair.first != -1)		// insert if present
-						LinetoAtomTest = LinetoAtom.insert( InsertionPair );
-					if (!LinetoAtomTest.second) {
-						// atom is already present in list, hence we have found a possible triangle
-						if ((Walker->nr < LinetoAtomTest.first->first) && (LinetoAtomTest.first->first < InsertionPair.first)) {	// check if numbering is in correct order for uniqueness, if so add
-							TriangleList[TriangleCount].x1 = Walker;
-							TriangleList[TriangleCount].x2 = LinetoAtomTest.first->second;
-							TriangleList[TriangleCount].x3 = InsertionPair.second;
-							TriangleCount++;
-						} 
-					}	//else {   // check if insertion was successful
-          // successfully inserted, hence nothing found yet
-        //} 
-				}
-			}
-		}
+    *out << Verbose(1) << "Examining line between " << (liner->first/AtomCount)+1 << " and " << (liner->first%AtomCount)+1 << "." << endl;
+    int enden[3][2];
+    enden[0][0] = (liner->first/AtomCount);
+    enden[0][1] = (liner->first%AtomCount);
+    if (enden[0][0] < enden[0][1]) {
+  		for (int endpoint=0;endpoint<2;endpoint++) {
+  			endpoints[0] = enden[0][endpoint];
+  			*out << Verbose(2) << "Current first endpoint is " << endpoints[0]+1 << "." << endl;
+  			
+  			LinesMap::iterator LineRangeStart = AllTriangleLines.lower_bound(endpoints[0]*AtomCount);
+  			for (LinesMap::iterator coach = LineRangeStart; (coach->first/AtomCount) == endpoints[0]; coach++) { // look at all line's other endpoint
+  		    enden[1][0] = (coach->first/AtomCount);
+  		    enden[1][1] = (coach->first%AtomCount);
+          //*out << Verbose(3) << "Line is " << coach->first << " with nr. " << coach->second << ": " << enden[1][0+1] << "<->" << enden[1][1]+1 << "." << endl;
+  			  endpoints[1] = (enden[1][0] != endpoints[0] ) ? enden[1][0] : enden[1][1];
+  			  if ((endpoints[1] > endpoints[0]) && (endpoints[1] != enden[0][(endpoint+1)%2])) {  // don't go back the very same line! 
+    			  *out << Verbose(3) << "Current second endpoint is " << endpoints[1]+1 << "." << endl;  
+    			  
+    	      LinesMap::iterator LineRangeStart2 = AllTriangleLines.lower_bound(endpoints[1]*AtomCount);
+    	      for (LinesMap::iterator coacher = LineRangeStart2; (coacher->first/AtomCount) == endpoints[1]; coacher++) { // look at all line's other endpoint
+    	        enden[2][0] = (coacher->first/AtomCount);
+    	        enden[2][1] = (coacher->first%AtomCount);
+              //*out << Verbose(3) << "Line is " << coacher->first << " with nr. " << coacher->second << ": " << enden[2][0]+1 << "<->" << enden[2][1]+1 << "." << endl;
+    	        endpoints[2] = (enden[2][0] != endpoints[1]) ? enden[2][0] : enden[2][1];
+    	        if ((endpoints[2] > endpoints[1]) && (endpoints[2] != endpoints[0])) {  // don't go back the very same line! 
+      	        *out << Verbose(4) << "Current third endpoint is " << endpoints[2]+1 << "." << endl;
+      	        
+      	        if (endpoints[2] == enden[0][(endpoint+1)%2]) { // jo, closing triangle
+                  *out << Verbose(3) << "MATCH: Triangle is ";
+                  for (int k=0;k<3;k++) {
+                    *out << endpoints[k]+1 << "\t";
+                    TriangleList[TriangleCount].x[k] = FindAtom(endpoints[k]); 
+                  }
+                  *out << endl;
+                  TriangleList[TriangleCount].nr = TriangleCount;
+                  TriangleCount++;
+      	        } else {
+      	          *out << Verbose(3) << "No match!" << endl;
+      	        }
+      	      }
+    	      }
+    			}
+  			}
+  		}
+    }
 	}
+  // listing created lines
+  *out << Verbose(2) << "The following triangles were created:";
+  for (int i=0;i<TriangleCount;i++) {
+    *out << " " << TriangleList[i].x[0]->Name << "<->" << TriangleList[i].x[1]->Name << "<->" << TriangleList[i].x[2]->Name; 
+  }
+  *out << endl;
 	*out << Verbose(2) << "I created " << TriangleCount << " triangles." << endl;
 
 	// 6. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+  *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
 	double volume = 0.;
 	double PyramidVolume = 0.;
@@ -1595,26 +1727,32 @@
 	double a,b,c;
 	for (int i=0; i<TriangleCount; i++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
-		x.CopyVector(&TriangleList[TriangleCount].x1->x);
-		x.SubtractVector(&TriangleList[TriangleCount].x2->x);
-		y.CopyVector(&TriangleList[TriangleCount].x1->x);
-		y.SubtractVector(&TriangleList[TriangleCount].x3->x);
-		a = sqrt(TriangleList[TriangleCount].x1->x.Distance(&TriangleList[TriangleCount].x2->x));
-		b = sqrt(TriangleList[TriangleCount].x1->x.Distance(&TriangleList[TriangleCount].x3->x));
-		c = sqrt(TriangleList[TriangleCount].x3->x.Distance(&TriangleList[TriangleCount].x2->x));
-		G =  c * (sin(x.Angle(&y))*a); // area of tesselated triangle
-		x.MakeNormalVector(&TriangleList[TriangleCount].x1->x, &TriangleList[TriangleCount].x2->x, &TriangleList[TriangleCount].x3->x);
-		y.CopyVector(&x);
-		x.Scale(CenterOfGravity->Projection(&x));
+		x.CopyVector(&TriangleList[i].x[0]->x);
+		x.SubtractVector(&TriangleList[i].x[1]->x);
+		y.CopyVector(&TriangleList[i].x[0]->x);
+		y.SubtractVector(&TriangleList[i].x[2]->x);
+		a = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[1]->x));
+		b = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[2]->x));
+		c = sqrt(TriangleList[i].x[2]->x.Distance(&TriangleList[i].x[1]->x));
+		G =  sqrt( ( (a*a+b*b+c*c)*(a*a+b*b+c*c) - 2*(a*a*a*a + b*b*b*b + c*c*c*c) )/16.); // area of tesselated triangle
+		x.MakeNormalVector(&TriangleList[i].x[0]->x, &TriangleList[i].x[1]->x, &TriangleList[i].x[2]->x);
+		x.Scale(TriangleList[i].x[1]->x.Projection(&x));
 		h = x.Norm(); // distance of CoG to triangle
 		PyramidVolume = (1./3.) * G * h;		// this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
+		*out << Verbose(2) << "Area of triangle is " << G << "A^2, height is " << h << " and the volume is " << PyramidVolume << "." << endl;
 		volume += PyramidVolume; 
 	}
 	*out << Verbose(1) << "The summed volume is " << volume << " " << (IsAngstroem ? "A^3" : "a_0^3") << "." << endl;
 
-	// free reference lists
+  // 7. translate all points back from CoG
+  *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
+	CenterOfGravity->Scale(-1);
+  Walker=start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    Walker->x.Translate(CenterOfGravity);
+  }
+
+  // free reference lists
 	delete[](AtomList);
-  for(int i=0;i<AtomCount;i++)
-    delete[](AtomLines[i]);
-  delete[](AtomLines);
 
 	return volume;
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision 233e33e9480b84b60d83d44696ee46721204ee12)
+++ src/molecules.hpp	(revision 49de64f43390d73584f072908bf52bbacde73c26)
@@ -40,13 +40,21 @@
 #define KeySet set<int>
 #define NumberValuePair pair<int, double>
-#define Graph map<KeySet, NumberValuePair, KeyCompare >
-#define GraphPair pair<KeySet, NumberValuePair >
+#define Graph map <KeySet, NumberValuePair, KeyCompare >
+#define GraphPair pair <KeySet, NumberValuePair >
 #define KeySetTestPair pair<KeySet::iterator, bool>
 #define GraphTestPair pair<Graph::iterator, bool>
 
-#define DistanceNrPair pair< double, atom* >
+#define DistanceNrPair pair < double, atom* >
+#define DistanceMap multimap < double, atom* >
+#define DistanceTestPair pair < DistanceMap::iterator, bool>
+
 #define Boundaries map <double, DistanceNrPair >
 #define BoundariesPair pair<double, DistanceNrPair >
-#define BoundariesTestPair pair<Boundaries::iterator, bool>
+#define BoundariesTestPair pair< Boundaries::iterator, bool>
+
+#define LinesMap map <int, int>
+#define LinesPair pair <int, int>
+#define LinesTestPair pair < LinesMap::iterator, bool>
+
 #define LinetoAtomMap map < int, atom * > 
 #define LinetoAtomPair pair < int, atom * > 
@@ -404,4 +412,5 @@
     char *configpath;
     char *configname;
+    bool FastParsing;
     
     private:
