Index: molecuilder/src/molecules.cpp
===================================================================
--- molecuilder/src/molecules.cpp	(revision 94b458e79385aa58ff62d8f35283066d1cb63942)
+++ molecuilder/src/molecules.cpp	(revision 32b6dc3848fc3a73112594cc2444b1e75c1a3357)
@@ -179,5 +179,4 @@
   bond *FirstBond = NULL, *SecondBond = NULL; // Other bonds in double bond case to determine "other" plane
   atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
-  int i;  // loop variable
   double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
   vector OrthoVector1, OrthoVector2;  // temporary vectors in coordination construction
@@ -253,5 +252,5 @@
     case 2:
       // determine two other bonds (warning if there are more than two other) plus valence sanity check
-      for (i=0;i<NumBond;i++) {
+      for (int i=0;i<NumBond;i++) {
         if (BondList[i] != TopBond) {
           if (FirstBond == NULL) {
@@ -312,5 +311,5 @@
       FirstOtherAtom->x.Zero();
       SecondOtherAtom->x.Zero();
-      for(i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
+      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
         FirstOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (sin(bondangle));
         SecondOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (-sin(bondangle));
@@ -319,5 +318,5 @@
       SecondOtherAtom->x.Scale(&BondRescale);
       //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
-      for(i=NDIM;i--;) { // and make relative to origin atom
+      for(int i=NDIM;i--;) { // and make relative to origin atom
         FirstOtherAtom->x.x[i] += TopOrigin->x.x[i]; 
         SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
@@ -439,5 +438,5 @@
   int NumberOfAtoms = 0; // atom number in xyz read
   int i, j; // loop variables
-  atom *first = NULL;  // pointer to added atom
+  atom *Walker = NULL;  // pointer to added atom
   char shorthand[3];  // shorthand for atom name
   ifstream xyzfile;   // xyz file
@@ -457,5 +456,5 @@
 
   for(i=0;i<NumberOfAtoms;i++){
-    first = new atom;
+    Walker = new atom;
     getline(xyzfile,line,'\n');
     istringstream *item = new istringstream(line);
@@ -466,12 +465,12 @@
     *item >> x[1];
     *item >> x[2];
-    first->type = elemente->FindElement(shorthand);
-    if (first->type == NULL) {
+    Walker->type = elemente->FindElement(shorthand);
+    if (Walker->type == NULL) {
       cerr << "Could not parse the element at line: '" << line << "', setting to H.";
-      first->type = elemente->FindElement(1);
+      Walker->type = elemente->FindElement(1);
     }
     for(j=NDIM;j--;)
-      first->x.x[j] = x[j];
-     AddAtom(first);  // add to molecule
+      Walker->x.x[j] = x[j];
+     AddAtom(Walker);  // add to molecule
     delete(item);
   }
@@ -710,17 +709,17 @@
 }; 
 
-/** Centers the center of gravity of the atoms at (0,0,0).
+/** Returns vector pointing to center of gravity.
  * \param *out output stream for debugging
- * \param *center return vector for translation vector
- */
-void molecule::CenterGravity(ofstream *out, vector *center)
-{
+ * \return pointer to center of gravity vector
+ */
+vector * molecule::DetermineCenterOfGravity(ofstream *out)
+{
+	atom *ptr = start->next;  // start at first in list
+	vector *a = new vector();
+	vector tmp;
   double Num = 0;
-  atom *ptr = start->next;  // start at first in list
-  vector tmp;
-  
-  for(int i=NDIM;i--;) // zero center vector
-    center->x[i] = 0.;
-    
+	
+	a->Zero();
+
   if (ptr != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
@@ -729,7 +728,25 @@
       tmp.CopyVector(&ptr->x);
       tmp.Scale(ptr->type->mass);  // scale by mass
-      center->AddVector(&tmp);       
-    }
-    center->Scale(-1./Num); // divide through total mass (and sign for direction)
+      a->AddVector(&tmp);       
+    }
+    a->Scale(-1./Num); // divide through total mass (and sign for direction)
+  }
+  *out << Verbose(1) << "Resulting center of gravity: ";
+  a->Output(out);
+  *out << endl;
+  return a;
+};
+
+/** Centers the center of gravity of the atoms at (0,0,0).
+ * \param *out output stream for debugging
+ * \param *center return vector for translation vector
+ */
+void molecule::CenterGravity(ofstream *out, vector *center)
+{
+  if (center == NULL) {
+    DetermineCenter(*center);
+    Translate(center);
+    delete(center);
+  } else {
     Translate(center);
   }
@@ -775,8 +792,8 @@
 };
 
-/** Determines center of gravity (yet not considering atom masses).
- * \param CenterOfGravity reference to return vector
- */
-void molecule::DetermineCenterOfGravity(vector &CenterOfGravity)
+/** Determines center of molecule (yet not considering atom masses).
+ * \param Center reference to return vector
+ */
+void molecule::DetermineCenter(vector &Center)
 {
   atom *Walker = start;
@@ -788,5 +805,5 @@
   
   do {
-    CenterOfGravity.Zero();
+    Center.Zero();
     flag = true;
     while (Walker->next != end) {
@@ -815,5 +832,5 @@
         TestVector.AddVector(&TranslationVector);
         TestVector.MatrixMultiplication(matrix);
-        CenterOfGravity.AddVector(&TestVector);
+        Center.AddVector(&TestVector);
         cout << Verbose(1) << "Vector is: ";
         TestVector.Output((ofstream *)&cout);
@@ -828,5 +845,5 @@
             TestVector.AddVector(&TranslationVector);
             TestVector.MatrixMultiplication(matrix);
-            CenterOfGravity.AddVector(&TestVector);
+            Center.AddVector(&TestVector);
             cout << Verbose(1) << "Hydrogen Vector is: ";
             TestVector.Output((ofstream *)&cout);
@@ -838,6 +855,110 @@
     }
   } while (!flag);
-  Free((void **)&matrix, "molecule::DetermineCenterOfGravity: *matrix");
-  CenterOfGravity.Scale(1./(double)AtomCount);
+  Free((void **)&matrix, "molecule::DetermineCenter: *matrix");
+  Center.Scale(1./(double)AtomCount);
+};
+
+/** Transforms/Rotates the given molecule into its principal axis system.
+ * \param *out output stream for debugging
+ * \param DoRotate whether to rotate (true) or only to determine the PAS.
+ */
+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;
+	    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.;
+    
+    // sum up inertia tensor
+    ptr = start;
+    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;
+	}
+	
+	// free everything
+	delete(CenterOfGravity);
+	gsl_vector_free(eval);
+	gsl_matrix_free(evec);
 };
 
@@ -1234,5 +1355,5 @@
   if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
     *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
-    Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
+    Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
     while (Subgraphs->next != NULL) {
       Subgraphs = Subgraphs->next;
@@ -1248,4 +1369,254 @@
   }
   return No;
+};
+
+/** Determines the volume of a cluster.
+ * Determines first the convex envelope, then tesselates it and calculates its volume.
+ * \param *out output stream for debugging
+ * \param IsAngstroem for the units on output
+ */
+double molecule::VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem) 
+{
+	atom *Walker = NULL;
+	
+	// 1. calculate center of gravity
+	vector *CenterOfGravity = DetermineCenterOfGravity(out);
+	
+	// 2. calculate distance of each atom and sort into hash table
+//	map<double, int> DistanceSet;
+//	
+//	Walker = start;
+//	while (Walker->next != end) {
+//		Walker = Walker->next;
+//		DistanceSet.insert( pair<double, int>(ptr->x.distance(CenterOfGravity), ptr->nr) ); 
+//	}
+	
+	// 3. Find all points on the boundary
+	Boundaries BoundaryPoints[NDIM];	// first is alpha, second is (r, nr)
+	BoundariesTestPair BoundaryTestPair;
+	vector AxisVector;
+	// 3a. Go through every axis
+	for (int axis=0; axis<NDIM; axis++)  {
+		AxisVector.Zero();
+		AxisVector.x[axis] = 1.;
+		// 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+		Walker = start;
+		while (Walker->next != end) {
+			Walker = Walker->next;
+			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) ) );
+			if (BoundaryTestPair.second) { // successfully inserted
+			} else { // same point exists, check first r, then distance of original vectors to centers of gravity
+				*out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
+				*out << Verbose(2) << "Present vector: ";
+				BoundaryTestPair.first->second.second->x.Output(out);
+				*out << endl;
+				*out << Verbose(2) << "New vector: ";
+				Walker->x.Output(out);
+				*out << endl;
+				double tmp = ProjectedVector.Norm(); 
+				if (tmp > BoundaryTestPair.first->second.first) {
+					BoundaryTestPair.first->second.first = tmp;
+					BoundaryTestPair.first->second.second = Walker; 
+					*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)) {
+						BoundaryTestPair.first->second.second = Walker;
+						*out << Verbose(2) << "Keeping new vector." << endl;
+					} else {
+						*out << Verbose(2) << "Keeping present vector." << endl;
+					}
+				} else {
+						*out << Verbose(2) << "Keeping present vector." << 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
+			flag = false;
+			Boundaries::iterator left = BoundaryPoints[axis].end();
+			Boundaries::iterator right = BoundaryPoints[axis].end();
+			for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+				// set neighbours correctly
+				if (runner == BoundaryPoints[axis].begin()) {
+					left = BoundaryPoints[axis].end();
+				}	else {
+					left = runner;
+				}
+        left--;
+        right = runner;
+        right++;
+				if (right == BoundaryPoints[axis].end()) {
+					right = BoundaryPoints[axis].begin();
+				}
+				// check distance
+				if (runner->second.first < (left->second.first + right->second.first)/2. ) {
+					// throw out point
+					BoundaryPoints[axis].erase(runner);
+					flag = true;
+				}
+			}
+		} while (flag);
+	}	
+	// 3d. put into boundary set only those points appearing in each of the NDIM sets
+	int *AtomList = new int[AtomCount];
+	for (int j=0; j<AtomCount; j++)	// reset list
+		AtomList[j] = 0;
+	for (int axis=0; axis<NDIM; axis++)  {	// increase whether it's present
+		for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+			AtomList[runner->second.second->nr]++;
+		}
+	}
+	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)
+				BoundaryPoints[axis].erase(runner);
+		}
+	}
+
+	*out << Verbose(2) << "I found " << BoundaryPointCount << " points on the convex boundary." << endl;
+	// now we have the whole set of edge points in the BoundaryList
+	
+	// 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;
+			}
+		}
+	}
+	
+	// then create each line
+  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++;
+  			}
+      }
+		}
+	}
+	*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;
+		int nr;
+	} TriangleList[2 * 3 * AtomCount];	// each line can be part in at most 2 triangles
+	int TriangleCount = 0;
+	for (int i=0; i<LineCount; i++) {	// 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(2) << "I created " << TriangleCount << " triangles." << endl;
+
+	// 6. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+	double volume = 0.;
+	double PyramidVolume = 0.;
+	double G,h;
+	vector x,y;
+	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));
+		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)
+		volume += PyramidVolume; 
+	}
+	*out << Verbose(1) << "The summed volume is " << volume << " " << (IsAngstroem ? "A^3" : "a_0^3") << "." << endl;
+
+	// free reference lists
+	delete[](AtomList);
+  for(int i=0;i<AtomCount;i++)
+    delete[](AtomLines[i]);
+  delete[](AtomLines);
+
+	return volume;
 };
 
@@ -1490,9 +1861,8 @@
  * We use the algorithm from [Even, Graph Algorithms, p.62].
  * \param *out output stream for debugging
- * \param ReturnStack true - return pointer to atom stack of separable components, false - return NULL
  * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found
  * \return list of each disconnected subgraph as an individual molecule class structure
  */
-MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize)
+MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, int *&MinimumRingSize)
 {
   class StackClass<atom *> *AtomStack;
@@ -2006,8 +2376,7 @@
  * \param *path path to file
  * \param *FragmentList empty, filled on return
- * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
  */
-bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList, bool IsAngstroem)
+bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
 {
   bool status = true;
@@ -2186,10 +2555,10 @@
 {
   ifstream File;
-  stringstream line;
+  stringstream filename;
   bool status = true;
   char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
   
-  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
-  File.open(line.str().c_str(), ios::out);
+  filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  File.open(filename.str().c_str(), ios::out);
   *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
   if (File != NULL) {
@@ -2493,10 +2862,10 @@
 
   // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
-  Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
+  Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
   // fill the bond structure of the individually stored subgraphs
   Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
 
   // ===== 3. if structure still valid, parse key set file and others =====
-  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList, configuration->GetIsAngstroem());
+  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
 
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
@@ -3421,5 +3790,4 @@
   KeyStack AtomStack;
   atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
-  KeySet::iterator runner;
   int RootKeyNr = FragmentSearch.Root->nr;
   int Counter = FragmentSearch.FragmentCounter;
@@ -3816,6 +4184,6 @@
 {
   Graph ***FragmentLowerOrdersList = NULL;
-  int Order, NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
-  int counter = 0;
+  int NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
+  int counter = 0, Order;
   int UpgradeCount = RootStack.size();
   KeyStack FragmentRootStack;
@@ -3831,6 +4199,4 @@
 
   // initialise the fragments structure
-  FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
-  FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
   FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
   FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
@@ -4059,6 +4425,6 @@
   if (result) {
     *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
-    DetermineCenterOfGravity(CenterOfGravity);
-    OtherMolecule->DetermineCenterOfGravity(OtherCenterOfGravity);
+    DetermineCenter(CenterOfGravity);
+    OtherMolecule->DetermineCenter(OtherCenterOfGravity);
     *out << Verbose(5) << "Center of Gravity: ";
     CenterOfGravity.Output(out);
