Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/Makefile.am	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -1,4 +1,4 @@
-SOURCE = atom.cpp bond.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
-HEADER = defs.hpp helpers.hpp molecules.hpp stackclass.hpp vector.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
+HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp stackclass.hpp vector.hpp
 
 bin_PROGRAMS = molecuilder joiner analyzer
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
+++ src/boundary.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -0,0 +1,939 @@
+#include "molecules.hpp"
+#include "boundary.hpp"
+
+// ======================================== Points on Boundary =================================
+
+BoundaryPointSet::BoundaryPointSet()
+{
+  LinesCount = 0;
+  Nr = -1;
+};
+
+BoundaryPointSet::BoundaryPointSet(atom *Walker)
+{
+  node = Walker;
+  LinesCount = 0;
+  Nr = Walker->nr;
+};
+
+BoundaryPointSet::~BoundaryPointSet()
+{
+  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  node = NULL;
+};
+
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line) 
+{
+  cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
+  if (line->endpoints[0] == this) {
+    lines.insert ( LinePair( line->endpoints[1]->Nr, line) );
+  } else {
+    lines.insert ( LinePair( line->endpoints[0]->Nr, line) );
+  }
+  LinesCount++;
+};
+
+ostream & operator << (ostream &ost, BoundaryPointSet &a) 
+{
+  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  return ost;
+};
+
+// ======================================== Lines on Boundary =================================
+
+BoundaryLineSet::BoundaryLineSet()
+{
+  for (int i=0;i<2;i++)
+    endpoints[i] = NULL;
+  TrianglesCount = 0;
+  Nr = -1;
+};
+
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+{
+  // set number
+  Nr = number;
+  // set endpoints in ascending order
+  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
+  // add this line to the hash maps of both endpoints
+  Point[0]->AddLine(this);
+  Point[1]->AddLine(this);
+  // clear triangles list
+  TrianglesCount = 0;
+  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+};
+
+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;
+      delete(endpoints[i]);
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
+};
+
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 
+{
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
+  triangles.insert ( TrianglePair( TrianglesCount, triangle) );
+  TrianglesCount++;
+};
+
+ostream & operator << (ostream &ost, BoundaryLineSet &a) 
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
+  return ost;
+};
+
+// ======================================== Triangles on Boundary =================================
+
+
+BoundaryTriangleSet::BoundaryTriangleSet()
+{
+  for (int i=0;i<3;i++) {
+    endpoints[i] = NULL;
+    lines[i] = NULL;
+  }
+  Nr = -1;
+};
+
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
+{
+  // set number
+  Nr = number;
+  // set lines
+  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
+  for (int i=0;i<3;i++) {
+    lines[i] = line[i];
+    lines[i]->AddTriangle(this);
+  }
+  // get ascending order of endpoints
+  map <int, class BoundaryPointSet * > OrderMap;
+  for(int i=0;i<3;i++)  // for all three lines
+    for (int j=0;j<2;j++) { // for both endpoints
+      OrderMap.insert ( pair <int, class BoundaryPointSet * >( line[i]->endpoints[j]->Nr, line[i]->endpoints[j]) );
+      // and we don't care whether insertion fails 
+    }
+  // set endpoints
+  int Counter = 0;
+  cout << Verbose(6) << " with end points ";
+  for (map <int, class BoundaryPointSet * >::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
+    endpoints[Counter] = runner->second;
+    cout << " " << *endpoints[Counter];
+    Counter++;
+  }
+  if (Counter < 3) {
+    cerr << "ERROR! We have a triangle with only two distinct endpoints!" << endl;
+    //exit(1);
+  }
+  cout << "." << endl;
+};
+
+BoundaryTriangleSet::~BoundaryTriangleSet()
+{
+  for (int i=0;i<3;i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    TriangleMap::iterator tester = lines[i]->triangles.begin();
+    tester++;
+    if (tester == lines[i]->triangles.end()) {
+      cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+      delete(lines[i]);
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+  }
+};
+
+void BoundaryTriangleSet::GetNormalVector(vector &NormalVector)
+{
+  // get normal vector
+  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x, &endpoints[2]->node->x);
+  
+  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+  if (endpoints[0]->node->x.Projection(&NormalVector) > 0)
+    NormalVector.Scale(-1.);
+};
+
+ostream & operator << (ostream &ost, BoundaryTriangleSet &a) 
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  return ost;
+};
+
+// ========================================== F U N C T I O N S =================================
+
+/** Finds the endpoint two lines are sharing.
+ * \param *line1 first line
+ * \param *line2 second line
+ * \return point which is shared or NULL if none
+ */
+class BoundaryPointSet * GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
+{
+  class BoundaryLineSet * lines[2] = {line1, line2};
+  class BoundaryPointSet *node = NULL;
+  map <int, class BoundaryPointSet * > OrderMap;
+  pair < map <int, class BoundaryPointSet * >::iterator, bool > OrderTest;
+  for(int i=0;i<2;i++)  // for both lines
+    for (int j=0;j<2;j++) { // for both endpoints
+      OrderTest = OrderMap.insert ( pair <int, class BoundaryPointSet * >( lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]) );
+      if (!OrderTest.second) { // if insertion fails, we have common endpoint
+        node = OrderTest.first->second;
+        cout << Verbose(5) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl;
+        j=2;
+        i=2;
+        break;
+      }
+    }
+  return node;
+};
+
+/** Determines the boundary points of a cluster.
+ * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
+ * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
+ * center and first and last point in the triple, it is thrown out.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure representing the cluster
+ */
+Boundaries * GetBoundaryPoints(ofstream *out, molecule *mol)
+{
+  atom *Walker = NULL;
+  PointMap PointsOnBoundary;
+  LineMap LinesOnBoundary;
+  TriangleMap TrianglesOnBoundary;
+  
+  *out << Verbose(1) << "Finding all boundary points." << endl;
+  Boundaries *BoundaryPoints = new Boundaries [NDIM];  // first is alpha, second is (r, nr)
+  BoundariesTestPair BoundaryTestPair;
+  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 = mol->start;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      vector ProjectedVector;
+      ProjectedVector.CopyVector(&Walker->x);
+      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 center 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.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;
+          } else {
+            *out << Verbose(2) << "Keeping present vector." << endl;
+          }
+        } else {
+            *out << Verbose(2) << "Keeping present vector." << endl;
+        }
+      }
+    }
+    // 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();
+      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
+        
+        // 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;
+          }
+        }
+      }
+    } while (flag);
+  } 
+  return BoundaryPoints;
+};
+
+/** Determines greatest diameters of a cluster defined by its convex envelope.
+ * Looks at lines parallel to one axis and where they intersect on the projected planes
+ * \param *out output stream for debugging
+ * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
+ * \param *mol molecule structure representing the cluster
+ * \param IsAngstroem whether we have angstroem or atomic units
+ * \return NDIM array of the diameters
+ */ 
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
+{
+  // get points on boundary of NULL was given as parameter
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
+  
+  // determine biggest "diameter" of cluster for each axis
+  Boundaries::iterator Neighbour, OtherNeighbour;
+  double *GreatestDiameter = new double[NDIM];
+  for(int i=0;i<NDIM;i++)
+    GreatestDiameter[i] = 0.;
+  double OldComponent, tmp, w1, w2;
+  vector DistanceVector, OtherVector;
+  int component, Othercomponent;
+  for(int axis=0;axis<NDIM;axis++) { // regard each projected plane
+    //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
+    for (int j=0;j<2;j++) { // and for both axis on the current plane
+      component = (axis+j+1)%NDIM;
+      Othercomponent = (axis+1+((j+1) & 1))%NDIM;
+      //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
+      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
+        // seek for the neighbours pair where the Othercomponent sign flips
+        Neighbour = runner;
+        Neighbour++;
+        if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
+          Neighbour = BoundaryPoints[axis].begin();
+        DistanceVector.CopyVector(&runner->second.second->x);
+        DistanceVector.SubtractVector(&Neighbour->second.second->x);
+        do {  // seek for neighbour pair where it flips
+          OldComponent = DistanceVector.x[Othercomponent]; 
+          Neighbour++;
+          if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
+            Neighbour = BoundaryPoints[axis].begin();
+          DistanceVector.CopyVector(&runner->second.second->x);
+          DistanceVector.SubtractVector(&Neighbour->second.second->x);
+          //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+        } while ((runner != Neighbour) && ( fabs( OldComponent/fabs(OldComponent) - DistanceVector.x[Othercomponent]/fabs(DistanceVector.x[Othercomponent]) ) < MYEPSILON)); // as long as sign does not flip
+        if (runner != Neighbour) {
+          OtherNeighbour = Neighbour;
+          if (OtherNeighbour == BoundaryPoints[axis].begin())  // make it wrap around
+            OtherNeighbour = BoundaryPoints[axis].end();
+          OtherNeighbour--;
+          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl; 
+          // now we have found the pair: Neighbour and OtherNeighbour
+          OtherVector.CopyVector(&runner->second.second->x);
+          OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
+          //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
+          //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
+          // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
+          w1 = fabs(OtherVector.x[Othercomponent]);
+          w2 = fabs(DistanceVector.x[Othercomponent]);
+          tmp = fabs((w1*DistanceVector.x[component] + w2*OtherVector.x[component])/(w1+w2));
+          // mark if it has greater diameter
+          //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
+          GreatestDiameter[component] = (GreatestDiameter[component] > tmp) ? GreatestDiameter[component] : tmp;
+        } //else
+          //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
+      }
+    }
+  }
+  *out << Verbose(0) << "RESULT: The biggest diameters are " << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and " << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "." << endl;
+
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[](BoundaryPoints);
+
+  return GreatestDiameter;
+};
+
+
+/** 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 *configuration needed for path to store convex envelope file
+ * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
+ * \param *mol molecule structure representing the cluster
+ */
+double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPtr, molecule *mol) 
+{
+  bool IsAngstroem = configuration->GetIsAngstroem();
+  atom *Walker = NULL;
+  struct Tesselation *TesselStruct = new Tesselation;
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G,h;
+  vector x,y;
+  double a,b,c;
+
+  // 1. calculate center of gravity
+  *out << endl;
+  vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
+  
+  // 2. translate all points into CoG
+  *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+  Walker = mol->start;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    Walker->x.Translate(CenterOfGravity);
+  }
+  
+  // 3. Find all points on the boundary
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
+  
+  // 4. fill the boundary point list
+  for (int axis=0;axis<NDIM;axis++)
+    for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+      TesselStruct->AddPoint(runner->second.second);
+    }
+
+  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
+  // now we have the whole set of edge points in the BoundaryList
+
+
+  // listing for debugging
+//  *out << Verbose(1) << "Listing PointsOnBoundary:";
+//  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+//    *out << " " << *runner->second;
+//  }
+//  *out << endl;
+  
+  // 5a. guess starting triangle
+  TesselStruct->GuessStartingTriangle(out);
+  
+  // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
+  TesselStruct->TesselateOnBoundary(out, configuration, mol);
+
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
+
+  // 6a. 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;
+  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
+    x.CopyVector(&runner->second->endpoints[0]->node->x);
+    x.SubtractVector(&runner->second->endpoints[1]->node->x);
+    y.CopyVector(&runner->second->endpoints[0]->node->x);
+    y.SubtractVector(&runner->second->endpoints[2]->node->x);
+    a = sqrt(runner->second->endpoints[0]->node->x.Distance(&runner->second->endpoints[1]->node->x));
+    b = sqrt(runner->second->endpoints[0]->node->x.Distance(&runner->second->endpoints[2]->node->x));
+    c = sqrt(runner->second->endpoints[2]->node->x.Distance(&runner->second->endpoints[1]->node->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(&runner->second->endpoints[0]->node->x, &runner->second->endpoints[1]->node->x, &runner->second->endpoints[2]->node->x);
+    x.Scale(runner->second->endpoints[1]->node->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 << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " << h << " and the volume is " << PyramidVolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+    volume += PyramidVolume; 
+  }
+  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+
+
+  // 7. translate all points back from CoG
+  *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
+  CenterOfGravity->Scale(-1);
+  Walker = mol->start;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    Walker->x.Translate(CenterOfGravity);
+  }
+
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[](BoundaryPoints);
+  
+  return volume;
+};
+
+
+/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
+ * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
+ * \param *out output stream for debugging
+ * \param *configuration needed for path to store convex envelope file
+ * \param *mol molecule structure representing the cluster
+ * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
+ * \param celldensity desired average density in final cell
+ */
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity)
+{
+  // transform to PAS
+  mol->PrincipalAxisSystem(out, true);
+  
+  // some preparations beforehand
+  bool IsAngstroem = configuration->GetIsAngstroem();
+  Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
+  double clustervolume;
+  if (ClusterVolume == 0)
+    clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
+  else 
+    clustervolume = ClusterVolume;
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
+  vector BoxLengths;
+  int repetition[NDIM] = {1, 1, 1};
+  int TotalNoClusters = 1;
+  for (int i=0;i<NDIM;i++)
+    TotalNoClusters *= repetition[i];
+
+  // sum up the atomic masses
+  double totalmass = 0.;
+  atom *Walker = mol->start;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    totalmass += Walker->type->mass;
+  }
+  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
+  
+  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass/clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
+  // solve cubic polynomial
+  *out << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
+  double cellvolume;
+  if (IsAngstroem)
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
+  else
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
+  *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
+  double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
+  *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  if (minimumvolume > cellvolume) {
+    cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
+    cout << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
+    for(int i=0;i<NDIM;i++)
+      BoxLengths.x[i] = GreatestDiameter[i];
+    mol->CenterEdge(out, &BoxLengths);
+  } else {
+    BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
+    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
+              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2] 
+              + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
+    BoxLengths.x[2] = minimumvolume - cellvolume;
+    double x0 = 0.,x1 = 0.,x2 = 0.;
+    if (gsl_poly_solve_cubic(BoxLengths.x[0],BoxLengths.x[1],BoxLengths.x[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
+      *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
+    else {
+      *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
+      x0 = x2;  // sorted in ascending order
+    }
+  
+    cellvolume = 1;
+    for(int i=0;i<NDIM;i++) {
+      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+      cellvolume *= BoxLengths.x[i];
+    }
+  
+    // set new box dimensions
+    *out << Verbose(0) << "Translating to box with these boundaries." << endl;
+    mol->CenterInBox((ofstream *)&cout, &BoxLengths);
+  }
+  // update Box of atoms by boundary
+  mol->SetBoxDimension(&BoxLengths);
+  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+};
+
+
+// =========================================================== class TESSELATION ===========================================
+
+/** Constructor of class Tesselation.
+ */
+Tesselation::Tesselation()
+{
+  PointsOnBoundaryCount = 0; 
+  LinesOnBoundaryCount = 0; 
+  TrianglesOnBoundaryCount = 0;
+};
+
+/** Constructor of class Tesselation.
+ * We have to free all points, lines and triangles.
+ */
+Tesselation::~Tesselation()
+{
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    delete(runner->second);
+  }
+};
+
+/** Gueses first starting triangle of the convex envelope.
+ * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
+ * \param *out output stream for debugging
+ * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
+ */ 
+void Tesselation::GuessStartingTriangle(ofstream *out)
+{
+  // 4b. create a starting triangle
+  // 4b1. create all distances
+  DistanceMultiMap DistanceMMap;
+  double distance;
+  for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+    for(PointMap::iterator sprinter = PointsOnBoundary.begin(); sprinter != PointsOnBoundary.end(); sprinter++) {
+      if (runner->first < sprinter->first) {
+        distance = runner->second->node->x.Distance(&sprinter->second->node->x);
+        DistanceMMap.insert( DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(runner,sprinter) ) );
+      }
+    }
+  }
+
+//    // listing distances
+//    *out << Verbose(1) << "Listing DistanceMMap:";
+//    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
+//      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
+//    }
+//    *out << endl;
+  
+  // 4b2. take three smallest distance that form a triangle
+  // we take the smallest distance as the base line
+  DistanceMultiMap::iterator baseline = DistanceMMap.begin();
+  BPS[0] = baseline->second.first->second;
+  BPS[1] = baseline->second.second->second;
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+
+  // take the second smallest as the second base line
+  DistanceMultiMap::iterator secondline = DistanceMMap.begin();
+  do {
+    secondline++;
+  } while (!(
+      (BPS[0] == secondline->second.first->second) && (BPS[1] != secondline->second.second->second) ||
+      (BPS[0] == secondline->second.second->second) && (BPS[1] != secondline->second.first->second) ||
+      (BPS[1] == secondline->second.first->second) && (BPS[0] != secondline->second.second->second) ||
+      (BPS[1] == secondline->second.second->second) && (BPS[0] != secondline->second.first->second)
+    ));
+  BPS[0] = secondline->second.first->second;
+  BPS[1] = secondline->second.second->second;
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+
+  // connection yields the third line (note: first and second endpoint are sorted!)
+  if (baseline->second.first->second == secondline->second.first->second) {
+    SetEndpointsOrdered(BPS, baseline->second.second->second, secondline->second.second->second);
+  } else if (baseline->second.first->second == secondline->second.second->second) {
+    SetEndpointsOrdered(BPS, baseline->second.second->second, secondline->second.first->second);
+  } else if (baseline->second.second->second == secondline->second.first->second) {
+    SetEndpointsOrdered(BPS, baseline->second.first->second, baseline->second.second->second);
+  } else if (baseline->second.second->second == secondline->second.second->second) {
+    SetEndpointsOrdered(BPS, baseline->second.first->second, baseline->second.first->second);
+  }
+  BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+  
+  // 4b3. insert created triangle
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+  for(int i=0;i<NDIM;i++) {
+    LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+    LinesOnBoundaryCount++;
+  }
+
+  *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+};
+
+
+/** Tesselates the convex envelope of a cluster from a single starting triangle.
+ * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
+ * 2 triangles. Hence, we go through all current lines:
+ * -# if the lines contains to only one triangle
+ * -# We search all points in the boundary
+ *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
+ *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
+ *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
+ *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
+ * \param *out output stream for debugging
+ * \param *configuration for IsAngstroem
+ * \param *mol the cluster as a molecule structure
+ */
+void Tesselation::TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol)
+{
+  bool flag;
+  PointMap::iterator winner;
+  class BoundaryPointSet *peak = NULL;
+  double SmallestAngle, TempAngle;
+  vector NormalVector, VirtualNormalVector, CenterVector, TempVector, PropagationVector;
+  LineMap::iterator LineChecker[2];
+  do {
+    flag = false;
+    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) 
+      if (baseline->second->TrianglesCount == 1) {
+        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl; 
+        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
+        SmallestAngle = M_PI;
+        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
+        // get peak point with respect to this base line's only triangle
+        for(int i=0;i<3;i++)
+          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
+            peak = BTS->endpoints[i];
+        *out << Verbose(3) << " and has peak " << *peak << "." << endl;
+        // normal vector of triangle
+        BTS->GetNormalVector(NormalVector);
+        *out << Verbose(4) << "NormalVector of base triangle is ";
+        NormalVector.Output(out);
+        *out << endl;
+        // offset to center of triangle
+        CenterVector.Zero();
+        for(int i=0;i<3;i++)
+          CenterVector.AddVector(&BTS->endpoints[i]->node->x);
+        CenterVector.Scale(1./3.);
+        *out << Verbose(4) << "CenterVector of base triangle is ";
+        CenterVector.Output(out);
+        *out << endl;
+        // vector in propagation direction (out of triangle)
+        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
+        TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
+        TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
+        PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
+        TempVector.CopyVector(&CenterVector);
+        TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x);  // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl; 
+        if (PropagationVector.Projection(&TempVector) > 0)  // make sure normal propagation vector points outward from baseline
+          PropagationVector.Scale(-1.);
+        *out << Verbose(4) << "PropagationVector of base triangle is ";
+        PropagationVector.Output(out);
+        *out << endl;
+        winner = PointsOnBoundary.end();
+        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) 
+          if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
+            *out << Verbose(3) << "Target point is " << *(target->second) << ":";
+            bool continueflag = true;
+            
+            VirtualNormalVector.CopyVector(&baseline->second->endpoints[0]->node->x);
+            VirtualNormalVector.AddVector(&baseline->second->endpoints[0]->node->x);
+            VirtualNormalVector.Scale(-1./2.);   // points now to center of base line
+            VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
+            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+            continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
+            if (!continueflag) {
+              *out << Verbose(4) << "Angle between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
+              continue;
+            } else
+              *out << Verbose(4) << "Angle between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
+            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
+            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
+  //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
+  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
+  //            else
+  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+  //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
+  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
+  //            else
+  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+            // check first endpoint (if any connecting line goes to target or at least not more than 1)
+            continueflag = continueflag && (( (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) || (LineChecker[0]->second->TrianglesCount == 1)));
+            if (!continueflag) {
+              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
+              continue;
+            }
+            // check second endpoint (if any connecting line goes to target or at least not more than 1)
+            continueflag = continueflag && (( (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) || (LineChecker[1]->second->TrianglesCount == 1)));
+            if (!continueflag) {
+              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
+              continue;
+            }
+            // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
+            continueflag = continueflag && (!(
+                ((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) 
+                && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak))
+               ));
+            if (!continueflag) {
+              *out << Verbose(4) << "Current target is peak!" << endl;
+              continue;
+            }
+            // in case NOT both were found
+            if (continueflag) {  // create virtually this triangle, get its normal vector, calculate angle
+              flag = true;
+              VirtualNormalVector.MakeNormalVector(&baseline->second->endpoints[0]->node->x, &baseline->second->endpoints[1]->node->x, &target->second->node->x);
+              // make it always point inward
+              if (baseline->second->endpoints[0]->node->x.Projection(&VirtualNormalVector) > 0)
+                VirtualNormalVector.Scale(-1.);
+              // calculate angle
+              TempAngle = NormalVector.Angle(&VirtualNormalVector);
+              *out << Verbose(4) << "NormalVector is ";
+              VirtualNormalVector.Output(out);
+              *out << " and the angle is " << TempAngle << "." << endl;
+              if (SmallestAngle > TempAngle) {  // set to new possible winner
+                SmallestAngle = TempAngle;
+                winner = target;
+              }
+            }
+          }
+        // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
+        if (winner != PointsOnBoundary.end()) {
+          *out << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;
+          // create the lins of not yet present
+          BLS[0] = baseline->second;
+          // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
+          LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
+          LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
+          if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
+            BPS[0] = baseline->second->endpoints[0];
+            BPS[1] = winner->second;
+            BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+            LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BLS[1]) );
+            LinesOnBoundaryCount++;
+          } else
+            BLS[1] = LineChecker[0]->second;
+          if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
+            BPS[0] = baseline->second->endpoints[1];
+            BPS[1] = winner->second;
+            BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+            LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BLS[2]) );
+            LinesOnBoundaryCount++;
+          } else
+            BLS[2] = LineChecker[1]->second;
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+          TrianglesOnBoundaryCount++;
+        } else {
+          *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
+        }
+  
+        // 5d. If the set of lines is not yet empty, go to 5. and continue
+      } else
+        *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
+  } while (flag);
+  
+  stringstream line;
+  line << configuration->configpath << "/" << CONVEXENVELOPE;
+  *out << Verbose(1) << "Storing convex envelope in tecplot data file " << line.str() << "." << endl;
+  ofstream output(line.str().c_str());
+  output << "TITLE = \"3D CONVEX SHELL\"" << endl;
+  output << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+  output << "ZONE T=\"TRIANGLES\", N=" << PointsOnBoundaryCount << ", E=" << TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+  int *LookupList = new int[mol->AtomCount];
+  for (int i=0;i<mol->AtomCount;i++)
+    LookupList[i] = -1;
+  
+  // print atom coordinates
+  *out << Verbose(2) << "The following triangles were created:";
+  int Counter = 1;
+  atom *Walker = NULL;
+  for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
+    Walker = target->second->node;
+    LookupList[Walker->nr] = Counter++;
+    output << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+  }
+  output << endl;
+    // print connectivity
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
+    output << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 
+  }
+  output.close();
+  delete[](LookupList);
+  *out << endl;
+};
+
+/** Adds an atom to the tesselation::PointsOnBoundary list.
+ * \param *Walker atom to add
+ */
+void Tesselation::AddPoint(atom *Walker)
+{
+  BPS[0] = new class BoundaryPointSet(Walker);
+  PointsOnBoundary.insert( PointPair(Walker->nr, BPS[0]) );
+  PointsOnBoundaryCount++;
+};
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
+++ src/boundary.hpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -0,0 +1,129 @@
+#ifndef BOUNDARY_HPP_
+#define BOUNDARY_HPP_
+
+class BoundaryPointSet;
+class BoundaryLineSet;
+class BoundaryTriangleSet;
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// STL headers
+#include <map>
+#include <set>
+#include <deque>
+
+#include <gsl/gsl_poly.h>
+
+#include "molecules.hpp"
+
+#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 PointMap map < int, class BoundaryPointSet * > 
+#define PointPair pair < int, class BoundaryPointSet * > 
+#define PointTestPair pair < PointMap::iterator, bool > 
+
+#define LineMap map < int, class BoundaryLineSet * > 
+#define LinePair pair < int, class BoundaryLineSet * > 
+#define LineTestPair pair < LinePair::iterator, bool > 
+
+#define TriangleMap map < int, class BoundaryTriangleSet * > 
+#define TrianglePair pair < int, class BoundaryTriangleSet * > 
+#define TriangleTestPair pair < TrianglePair::iterator, bool > 
+
+#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
+#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
+
+template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
+{
+  if (endpoint1->Nr < endpoint2->Nr) {
+    endpoints[0] = endpoint1;
+    endpoints[1] = endpoint2;
+  } else {
+    endpoints[0] = endpoint2;
+    endpoints[1] = endpoint1;
+  }
+};
+
+class BoundaryPointSet {
+  public:
+    BoundaryPointSet();
+    BoundaryPointSet(atom *Walker);
+    ~BoundaryPointSet();
+    
+    void AddLine(class BoundaryLineSet *line);
+    
+    LineMap lines;
+    int LinesCount;
+    atom *node;
+    int Nr;
+};
+
+class BoundaryLineSet {
+  public:
+    BoundaryLineSet();
+    BoundaryLineSet(class BoundaryPointSet *Point[2], int number);
+    ~BoundaryLineSet();
+
+    void AddTriangle(class BoundaryTriangleSet *triangle);
+
+    class BoundaryPointSet *endpoints[2];
+    TriangleMap triangles;
+    int TrianglesCount;
+    int Nr;
+};
+
+class BoundaryTriangleSet {
+  public:
+    BoundaryTriangleSet();
+    BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
+    ~BoundaryTriangleSet();
+    
+    void GetNormalVector(vector &NormalVector);
+    
+    class BoundaryPointSet *endpoints[3];
+    class BoundaryLineSet *lines[3];
+    int Nr;
+};
+
+class Tesselation {
+  public:
+    
+    Tesselation();
+    ~Tesselation();
+    
+    void TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol);
+    void GuessStartingTriangle(ofstream *out);
+    void AddPoint(atom * Walker);
+    
+    PointMap PointsOnBoundary;
+    LineMap LinesOnBoundary;
+    TriangleMap TrianglesOnBoundary;
+    class BoundaryPointSet *BPS[2];
+    class BoundaryLineSet *BLS[3];
+    class BoundaryTriangleSet *BTS;
+    int PointsOnBoundaryCount;
+    int LinesOnBoundaryCount;
+    int TrianglesOnBoundaryCount;
+};
+
+
+ostream & operator << (ostream &ost, BoundaryPointSet &a);
+ostream & operator << (ostream &ost, BoundaryLineSet &a);
+ostream & operator << (ostream &ost, BoundaryTriangleSet &a);
+
+
+double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPoints, molecule *mol);
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
+
+
+#endif /*BOUNDARY_HPP_*/
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/builder.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -52,4 +52,5 @@
 #include "helpers.hpp"
 #include "molecules.hpp"
+#include "boundary.hpp"
 
 /********************************************** Submenu routine **************************************/
@@ -484,5 +485,5 @@
  * \param *mol the molecule with all the atoms
  */
-static void MeasureAtoms(periodentafel *periode, molecule *mol)
+static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
 {
   atom *first, *second, *third;
@@ -496,4 +497,6 @@
   cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
   cout << Verbose(0) << " c - calculate bond angle" << endl;
+  cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
+  cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
   cout << Verbose(0) << "all else - go back" << endl;
   cout << Verbose(0) << "===============================================" << endl;
@@ -553,4 +556,17 @@
       cout << Verbose(0) << (acos(x.ScalarProduct((const vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;         
       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;
+    case 'e':
+    	cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+    	VolumeOfConvexEnvelope((ofstream *)&cout, configuration, NULL, mol);
+    	break;
   }
 };
@@ -720,4 +736,5 @@
   int ExitFlag = 0;
   int j;
+  double volume = 0.;
   enum ConfigStatus config_present = absent;
   clock_t start,end;
@@ -742,11 +759,16 @@
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
+            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
             cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
+            cout << "\t-m\tAlign in PAS with greatest EV along z axis." << endl;
+            cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
+            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
+            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 
             cout << "\t-v/-V\t\tGives version information." << endl;
             cout << "Note: config files must not begin with '-' !" << endl;
@@ -768,4 +790,7 @@
             argptr+=1;
             break;
+          case 'n':
+            cout << "I won't parse trajectories." << endl;
+            configuration.FastParsing = true;
           default:   // no match? Step on
             argptr++;
@@ -780,6 +805,8 @@
       cout << Verbose(0) << "Element list loaded successfully." << endl;
       periode->Output((ofstream *)&cout);
-    } else
+    } else {
       cout << Verbose(0) << "Element list loading failed." << endl;
+      return 1;
+    }
     
     // 3. Find config file name and parse if possible
@@ -945,6 +972,99 @@
               }
               break;
+            case 'm':
+              ExitFlag = 1;
+              j = atoi(argv[argptr++]);
+              if ((j<0) || (j>1)) {
+                cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
+                j = 0;
+              }
+              if (j)
+                cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
+              else
+                cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+              mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
+              break;
+            case 'o':
+              ExitFlag = 1;
+              cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+              VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
+              break;
+            case 'U':
+              volume = atof(argv[argptr++]);
+              cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;  
+            case 'u':
+              {
+                double density;
+                ExitFlag = 1;
+                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
+                density = atof(argv[argptr++]);
+                if (density < 1.0) {
+                  cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
+                  density = 1.3;
+                }
+//                for(int i=0;i<NDIM;i++) {
+//                  repetition[i] = atoi(argv[argptr++]);
+//                  if (repetition[i] < 1)
+//                    cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
+//                  repetition[i] = 1;
+//                }
+                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);
+              }
+              break;
+            case 'd':
+              for (int axis = 1; axis <= NDIM; axis++) {
+                int faktor = atoi(argv[argptr++]);
+                int count;
+                element ** Elements;
+                vector ** Vectors;
+                if (faktor < 1) {
+                  cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
+                  faktor = 1;
+                }
+                mol->CountAtoms((ofstream *)&cout);  // recount atoms
+                if (mol->AtomCount != 0) {  // if there is more than none
+                  count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
+                  Elements = new element *[count];
+                  Vectors = new vector *[count];
+                  j = 0;
+                  first = mol->start;
+                  while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
+                    first = first->next;
+                    Elements[j] = first->type;
+                    Vectors[j] = &first->x;
+                    j++;
+                  }
+                  if (count != j)
+                    cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+                  x.Zero();
+                  y.Zero();
+                  y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+                  for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
+                    x.AddVector(&y); // per factor one cell width further
+                    for (int k=count;k--;) { // go through every atom of the original cell
+                      first = new atom(); // create a new body
+                      first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
+                      first->x.AddVector(&x);      // translate the coordinates
+                      first->type = Elements[k];  // insert original element
+                      mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
+                    }
+                  }
+                  // free memory
+                  delete[](Elements);
+                  delete[](Vectors);
+                  // correct cell size
+                  if (axis < 0) { // if sign was negative, we have to translate everything
+                    x.Zero();
+                    x.AddVector(&y);
+                    x.Scale(-(faktor-1));
+                    mol->Translate(&x);
+                  }
+                  mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+                }
+              }
+              break;
             default:   // no match? Step on
-              argptr++;
+              if (argv[argptr][0] != '-') // if it started with a '-' we've already made a step!
+                argptr++;
               break;
           }
@@ -1151,5 +1271,5 @@
 
       case 'l': // measure distances or angles
-        MeasureAtoms(periode, mol);
+        MeasureAtoms(periode, mol, &configuration);
         break;
 
@@ -1164,5 +1284,5 @@
         mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
         //mol->CreateListOfBondsPerAtom((ofstream *)&cout);
-        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
+        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, MinimumRingSize);
         while (Subgraphs->next != NULL) {
           Subgraphs = Subgraphs->next;
@@ -1220,5 +1340,5 @@
   
   // save element data base
-  if (periode->StorePeriodentafel()) //ElementsFileName
+  if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
     cout << Verbose(0) << "Saving of elements.db successful." << endl;
   else
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/config.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -24,4 +24,5 @@
   configname[0]='\0';
   
+  FastParsing = false;
   ProcPEGamma=8;
   ProcPEPsi=1;
@@ -455,8 +456,8 @@
   string zeile;
   string dummy;
-  element *elementhash[128];
-  char name[128];
-  char keyword[128];
-  int Z, No;
+  element *elementhash[MAX_ELEMENTS];
+  char name[MAX_ELEMENTS];
+  char keyword[MAX_ELEMENTS];
+  int Z, No[MAX_ELEMENTS];
   int verbose = 0;
   
@@ -632,20 +633,28 @@
   if (!ParseForParameter(verbose,file,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
     config::StructOpt = 0;
+  // prescan number of ions per type
+  cout << Verbose(0) << "Prescanning ions per type: " << endl;
   for (int i=0; i < config::MaxTypes; i++) {
     sprintf(name,"Ion_Type%i",i+1);
-    ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No, 1, critical);
+    ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
     ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
     elementhash[i] = periode->FindElement(Z);
-    for(int j=0;j<No;j++) {
+    cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl; 
+  }
+  for (int i=0; i < config::MaxTypes; i++) {
+    sprintf(name,"Ion_Type%i",i+1);
+    for(int j=0;j<No[i];j++) {
       int repetition = 1; // which repeated keyword shall be read
       sprintf(keyword,"%s_%i",name, j+1);
       atom *neues = new atom();
       // then parse for each atom the coordinates as often as present
-      while (ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition, (repetition == 1) ? critical : optional) &&
-        ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition, (repetition == 1) ? critical : optional) &&
-        ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition, (repetition == 1) ? critical : optional))
-        repetition++;
-      repetition--;
-      //cout << "Found " << repetition << " times of the keyword " << keyword << "." << endl;
+      if (!FastParsing) {
+        while (ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition, (repetition == 1) ? critical : optional) &&
+          ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition, (repetition == 1) ? critical : optional) &&
+          ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition, (repetition == 1) ? critical : optional))
+          repetition++;
+        repetition--;
+        //cout << "Found " << repetition << " times of the keyword " << keyword << "." << endl;
+      }
       ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
       ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
@@ -656,5 +665,5 @@
       if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
         neues->v.x[1] = 0.;
-      if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
+      if(!ParseForParameter(verbose,file, keyword, (int)FastParsing, 7, 1, double_type, &neues->v.x[2], repetition,optional))
         neues->v.x[2] = 0.;
       neues->type = elementhash[i]; // find element type
Index: src/defs.hpp
===================================================================
--- src/defs.hpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/defs.hpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -43,18 +43,24 @@
 
 // various standard filenames
-#define DEFAULTCONFIG "main_pcp_linux"
-#define KEYSETFILE "KeySets.dat"
-#define ADJACENCYFILE "Adjacency.dat"
-#define TEFACTORSFILE "TE-Factors.dat"
-#define FORCESFILE "Forces-Factors.dat"
-#define ORDERATSITEFILE "OrderAtSite.dat"
-#define ENERGYPERFRAGMENT "EnergyPerFragment"
-#define FRAGMENTPREFIX "BondFragment"
-#define STANDARDCONFIG "unknown.conf"
-#define STANDARDELEMENTSDB "elements.db"
-#define STANDARDVALENCEDB "valence.db"
-#define STANDARDORBITALDB "orbitals.db"
-#define STANDARDHBONDDISTANCEDB "Hbonddistance.db"
-#define STANDARDHBONDANGLEDB "Hbondangle.db"
+#define DEFAULTCONFIG "main_pcp_linux"    //!< default filename of config file
+#define CONVEXENVELOPE "ConvexEnvelope.dat"    //!< default filename of convex envelope tecplot data file
+#define KEYSETFILE "KeySets.dat"    //!< default filename of BOSSANOVA key sets file
+#define ADJACENCYFILE "Adjacency.dat"    //!< default filename of BOSSANOVA adjacancy file
+#define TEFACTORSFILE "TE-Factors.dat"    //!< default filename of BOSSANOVA total energy factors file
+#define FORCESFILE "Forces-Factors.dat"    //!< default filename of BOSSANOVA force factors file
+#define ORDERATSITEFILE "OrderAtSite.dat"    //!< default filename of BOSSANOVA Bond Order at each atom file
+#define ENERGYPERFRAGMENT "EnergyPerFragment"    //!< default filename of BOSSANOVA Energy contribution Per Fragment file
+#define FRAGMENTPREFIX "BondFragment"    //!< default filename prefix of BOSSANOVA fragment config and directories 
+#define STANDARDCONFIG "unknown.conf"    //!< default filename of standard config file
+#define STANDARDELEMENTSDB "elements.db"    //!< default filename of elements data base with masses, Z, VanDerWaals radii, ...
+#define STANDARDVALENCEDB "valence.db"    //!< default filename of valence number per element database
+#define STANDARDORBITALDB "orbitals.db"    //!< default filename of orbitals per element database
+#define STANDARDHBONDDISTANCEDB "Hbonddistance.db"    //!< default filename of typial bond distance to hydrogen database
+#define STANDARDHBONDANGLEDB "Hbondangle.db"    //!< default filename of typial bond angle to hydrogen database
+
+// some values
+#define SOLVENTDENSITY_A 0.6022142
+#define SOLVENTDENSITY_a0 0.089238936
+
 
 #define UPDATECOUNT 10  //!< update ten sites per BOSSANOVA interval
Index: src/helpers.hpp
===================================================================
--- src/helpers.hpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/helpers.hpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -118,6 +118,4 @@
 };
 
-
-
 /******************************** Some templates for list management ***********************************/
 
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/moleculelist.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -240,8 +240,8 @@
     *out << endl;
     // prepare output of config file
-    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.conf", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    strcpy(PathBackup, configuration->configpath);
-    sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    //strcpy(PathBackup, configuration->configpath);
+    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     
     // center on edge
@@ -251,5 +251,5 @@
     for (int k=0;k<NDIM;k++) {
       j += k+1;
-      BoxDimension.x[k] = 5.;
+      BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
       ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
     }
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/molecules.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -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;
@@ -1249,5 +1370,4 @@
   return No;
 };
-
 /** Returns Shading as a char string.
  * \param color the Shading
@@ -1453,5 +1573,5 @@
             for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
               NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-	          //*out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+	          *out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
 	          if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check NoBonds of other atom
 	            // count valence of second partner
@@ -1459,5 +1579,5 @@
               for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
                 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
-              //*out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+              *out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
 	            if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) // increase bond degree by one
                 ListOfBondsPerAtom[Walker->nr][i]->BondDegree++;
@@ -1471,12 +1591,12 @@
 	  *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << "." << 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 << 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 << endl;
   } else
   	*out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl; 
@@ -1490,9 +1610,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 +2125,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 +2304,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 +2611,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)
@@ -3013,5 +3131,5 @@
 //          *out << ", who has no son in this fragment molecule." << endl; 
 #ifdef ADDHYDROGEN
-//          *out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
+          //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
           Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem);
 #endif
@@ -3421,5 +3539,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 +3933,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 +3948,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 +4174,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);
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/molecules.hpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -13,4 +13,5 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_matrix.h>
+#include <gsl/gsl_eigen.h>
 #include <gsl/gsl_heapsort.h>
 
@@ -39,6 +40,6 @@
 #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>
@@ -265,14 +266,17 @@
   void CenterEdge(ofstream *out, vector *max); 
   void CenterOrigin(ofstream *out, vector *max); 
-  void CenterGravity(ofstream *out, vector *max); 
+  void CenterGravity(ofstream *out, vector *max);
   void Translate(const vector *x);
   void Mirror(const vector *x);
   void Align(vector *n);
   void Scale(double **factor);
-  void DetermineCenterOfGravity(vector &CenterOfGravity);
+  void DetermineCenter(vector &center);
+  vector * DetermineCenterOfGravity(ofstream *out);
   void SetBoxDimension(vector *dim);
   double * ReturnFullMatrixforSymmetric(double *cell_size);
   void ScanForPeriodicCorrection(ofstream *out);
-
+	void PrincipalAxisSystem(ofstream *out, bool DoRotate);
+	double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
+	
   bool CheckBounds(const vector *x) const;
   void GetAlignVector(struct lsq_params * par) const;
@@ -283,5 +287,5 @@
   
   // Graph analysis
-  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize);
+  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, int *&MinimumRingSize);
   void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
   bond * FindNextUnused(atom *vertex);
@@ -303,5 +307,5 @@
   bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
   bool StoreOrderAtSiteFile(ofstream *out, char *path);
-  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList, bool IsAngstroem);
+  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
   bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
   bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
@@ -392,4 +396,5 @@
     char *configpath;
     char *configname;
+    bool FastParsing;
     
     private:
Index: src/periodentafel.cpp
===================================================================
--- src/periodentafel.cpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/periodentafel.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -217,6 +217,7 @@
 
   // fill valence DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -226,5 +227,5 @@
     	infile >> FindElement((int)tmp)->Valence;
     	infile >> ws;
-    	//cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->Valence << " valence electrons." << endl;
+    	//cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->Valence << " valence electrons." << endl;
     }
     infile.close();
@@ -234,6 +235,7 @@
 
   // fill valence DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -243,5 +245,5 @@
       infile >> FindElement((int)tmp)->NoValenceOrbitals;
       infile >> ws;
-      //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
     }
     infile.close();
@@ -251,6 +253,7 @@
   
   // fill H-BondDistance DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -263,5 +266,5 @@
       infile >> ptr->HBondDistance[2];
     	infile >> ws;
-      //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
     }
     infile.close();
@@ -271,6 +274,7 @@
   
   // fill H-BondAngle DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -289,4 +293,7 @@
     otherstatus = false;
   
+  if (!otherstatus)
+    cerr << "ERROR: Something went wrong while parsing the databases!" << endl;
+  
   return status;
 };
@@ -298,7 +305,8 @@
   bool result = true;
   ofstream f;
-  
-  if (filename == NULL)
-    f.open("elements.db");
+  char file[MAXSTRINGSIZE];
+  
+  if (filename == STANDARDELEMENTSDB)
+    f.open(file);
   else
     f.open(filename);
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/vector.cpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -13,4 +13,8 @@
  */
 vector::vector() { x[0] = x[1] = x[2] = 0.; };
+
+/** Constructor of class vector.
+ */
+vector::vector(double x1, double x2, double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
 
 /** Desctructor of class vector.
@@ -111,4 +115,16 @@
 };
 
+/** projects this vector onto plane defined by \a *y.
+ * \param *y array to normal vector of plane
+ * \return \f$\langle x, y \rangle\f$
+ */
+void vector::ProjectOntoPlane(const vector *y)
+{
+  vector tmp;
+  tmp.CopyVector(y);
+  tmp.Scale(Projection(y));
+  this->SubtractVector(&tmp);
+};
+
 /** Calculates the projection of a vector onto another \a *y.
  * \param *y array to second vector
@@ -117,5 +133,5 @@
 double vector::Projection(const vector *y) const
 {
-  return (ScalarProduct(y)/Norm()/y->Norm());
+  return (ScalarProduct(y));
 };
 
@@ -150,11 +166,28 @@
 };
 
+/** Zeros all components of this vector.
+ */
+void vector::One(double one)
+{
+  for (int i=NDIM;i--;)
+    this->x[i] = one;
+};
+
+/** Initialises all components of this vector.
+ */
+void vector::Init(double x1, double x2, double x3)
+{ 
+  x[0] = x1;
+  x[1] = x2;
+  x[2] = x3;
+};
+
 /** Calculates the angle between this and another vector.
  * \param *y array to second vector
- * \return \f$\frac{\langle x, y \rangle}{|x||y|}\f$
+ * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
  */
 double vector::Angle(vector *y) const
 {
-  return (this->ScalarProduct(y)/(this->Norm()*y->Norm()));
+  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
 };
 
@@ -354,5 +387,5 @@
 {
   double projection;
-  projection = ScalarProduct(n)/((vector *)n)->ScalarProduct(n);    // remove constancy from n (keep as logical one)
+  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
   // withdraw projected vector twice from original one
   cout << Verbose(1) << "Vector: ";
@@ -385,10 +418,10 @@
     return false;
   }
-  cout << Verbose(4) << "relative, first plane coordinates:";
-  x1.Output((ofstream *)&cout);
-  cout << endl;
-  cout << Verbose(4) << "second plane coordinates:";
-  x2.Output((ofstream *)&cout);
-  cout << endl;
+//  cout << Verbose(4) << "relative, first plane coordinates:";
+//  x1.Output((ofstream *)&cout);
+//  cout << endl;
+//  cout << Verbose(4) << "second plane coordinates:";
+//  x2.Output((ofstream *)&cout);
+//  cout << endl;
 
   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
@@ -419,10 +452,10 @@
     return false;
   }
-  cout << Verbose(4) << "relative, first plane coordinates:";
-  x1.Output((ofstream *)&cout);
-  cout << endl;
-  cout << Verbose(4) << "second plane coordinates:";
-  x2.Output((ofstream *)&cout);
-  cout << endl;
+//  cout << Verbose(4) << "relative, first plane coordinates:";
+//  x1.Output((ofstream *)&cout);
+//  cout << endl;
+//  cout << Verbose(4) << "second plane coordinates:";
+//  x2.Output((ofstream *)&cout);
+//  cout << endl;
 
   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
@@ -458,5 +491,5 @@
  * \return true - success, false - failure (null vector given)
  */
-bool vector::GetOneNormalVector(const vector *vector)
+bool vector::GetOneNormalVector(const vector *GivenVector)
 {
   int Components[NDIM]; // contains indices of non-zero components
@@ -466,5 +499,5 @@
 
   cout << Verbose(4);
-  vector->Output((ofstream *)&cout);
+  GivenVector->Output((ofstream *)&cout);
   cout << endl;
   for (j=NDIM;j--;)
@@ -472,5 +505,5 @@
   // find two components != 0
   for (j=0;j<NDIM;j++)
-    if (fabs(vector->x[j]) > MYEPSILON)
+    if (fabs(GivenVector->x[j]) > MYEPSILON)
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
@@ -479,9 +512,9 @@
     case 3:  // threecomponent system
     case 2:  // two component system
-      norm = sqrt(1./(vector->x[Components[1]]*vector->x[Components[1]]) + 1./(vector->x[Components[0]]*vector->x[Components[0]]));
+      norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
       x[Components[2]] = 0.;
       // in skp both remaining parts shall become zero but with opposite sign and third is zero
-      x[Components[1]] = -1./vector->x[Components[1]] / norm;
-      x[Components[0]] = 1./vector->x[Components[0]] / norm;
+      x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
+      x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
       return true;
       break;
@@ -498,4 +531,20 @@
 };
 
+/** Determines paramter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
+ * \param *A first plane vector
+ * \param *B second plane vector
+ * \param *C third plane vector
+ * \return scaling parameter for this vector
+ */
+double vector::CutsPlaneAt(vector *A, vector *B, vector *C)
+{
+//  cout << Verbose(3) << "For comparison: ";
+//  cout << "A " << A->Projection(this) << "\t"; 
+//  cout << "B " << B->Projection(this) << "\t"; 
+//  cout << "C " << C->Projection(this) << "\t"; 
+//  cout << endl;
+  return A->Projection(this);
+};
+
 /** Creates a new vector as the one with least square distance to a given set of \a vectors.
  * \param *vectors set of vectors
@@ -519,5 +568,5 @@
      gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *s = NULL;
-   gsl_vector *ss, *x;
+   gsl_vector *ss, *y;
    gsl_multimin_function minex_func;
  
@@ -528,5 +577,5 @@
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
-   x = gsl_vector_alloc (np);
+   y = gsl_vector_alloc (np);
  
    /* Set all step sizes to 1 */
@@ -538,5 +587,5 @@
 	
 	 for (i=NDIM;i--;)
-		gsl_vector_set(x, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
+		gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
 	 
    /* Initialize method and iterate */
@@ -546,5 +595,5 @@
  
    s = gsl_multimin_fminimizer_alloc (T, np);
-   gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
+   gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
  
    do
@@ -575,5 +624,5 @@
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
-   gsl_vector_free(x);
+   gsl_vector_free(y);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free (s);
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision d250b29e2b8e1c054b7f93c4f475f36dfb7a3a33)
+++ src/vector.hpp	(revision 8a3e457fdc8672973ca81864d2f77e092e907b83)
@@ -10,4 +10,5 @@
 
   vector();
+  vector(double x1, double x2, double x3);
   ~vector();
 
@@ -23,5 +24,8 @@
   void CopyVector(const vector *y);
   void RotateVector(const vector *y, const double alpha);
+  void ProjectOntoPlane(const vector *y);
   void Zero(); 
+  void One(double one);
+  void Init(double x1, double x2, double x3);
   void Normalize();
   void Translate(const vector *x);
@@ -35,4 +39,5 @@
   void LinearCombinationOfVectors(const vector *x1, const vector *x2, const vector *x3, double *factors);
   
+  double CutsPlaneAt(vector *A, vector *B, vector *C);
   bool GetOneNormalVector(const vector *x1);
   bool MakeNormalVector(const vector *y1);
