Changeset 357fba for src

Timestamp:

Aug 3, 2009, 2:48:42 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

a80fbdf, caf4ba

Parents:

2319ed

Message:

Huge refactoring of Tesselation routines, but not finished yet.

• new file tesselation.cpp with all of classes tesselation, Boundary..Set and CandidatesForTesselationOB
• new file tesselationhelper.cpp with all auxiliary functions.
• boundary.cpp just contains super functions, combininb molecule and Tesselation pointers
• new pointer molecule::TesselStruct
• PointMap, LineMap, TriangleMap DistanceMap have been moved from molecules.hpp to tesselation.hpp
• new abstract class PointCloud and TesselPoint
• atom inherits TesselPoint
• molecule inherits PointCloud (i.e. a set of TesselPoints) and implements all virtual functions for the chained list
• TriangleFilesWritten is thrown out, intermediate steps are written in find_nonconvex_border and not in find_next_triangle()
• LinkedCell class uses TesselPoint as its nodes, i.e. as long as any class inherits TesselPoint, it may make use of LinkedCell as well and a PointCloud is used to initialize
• class atom and bond definitions have been moved to own header files

NOTE: This is not bugfree yet. Tesselation of heptan produces way too many triangles, but runs without faults or leaks.

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• atom.cpp (modified) (5 diffs)
• atom.hpp (added)
• bond.cpp (modified) (1 diff)
• bond.hpp (added)
• boundary.cpp (modified) (11 diffs)
• boundary.hpp (modified) (3 diffs)
• builder.cpp (modified) (1 diff)
• config.cpp (modified) (9 diffs)
• ellipsoid.cpp (modified) (5 diffs)
• ellipsoid.hpp (modified) (1 diff)
• linkedcell.cpp (modified) (7 diffs)
• linkedcell.hpp (modified) (3 diffs)
• molecules.cpp (modified) (2 diffs)
• molecules.hpp (modified) (8 diffs)
• tesselation.cpp (added)
• tesselation.hpp (added)
• tesselationhelpers.cpp (added)
• tesselationhelpers.hpp (added)
• vector.hpp (modified) (1 diff)

src/Makefile.am

@@ -1 @@
-SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp vector.cpp verbose.cpp
-HEADER = boundary.hpp defs.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp vector.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
+HEADER = atom.hpp bond.hpp boundary.hpp defs.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp
 
 bin_PROGRAMS = molecuilder joiner analyzer

src/atom.cpp

@@ -5 +5 @@
  */
 
-#include "molecules.hpp"
+#include "atom.hpp"
 
 /************************************* Functions for class atom *************************************/
@@ -13 +13 @@
 atom::atom()
 {
-  Name = NULL;
+  father = this;  // generally, father is itself
   previous = NULL;
   next = NULL;
-  father = this;  // generally, father is itself
   Ancestor = NULL;
   type = NULL;
   sort = NULL;
   FixedIon = 0;
-  nr = -1;
   GraphNr = -1;
   ComponentNr = NULL;
@@ -29 +27 @@
   AdaptiveOrder = 0;
   MaxOrder = false;
+  // set LCNode::Vector to our Vector
+  node = &x;
 };
 
 /** Constructor of class atom.
- *
  */
 atom::atom(atom *pointer)
@@ -41 +40 @@
   father = this;  // generally, father is itself
   Ancestor = NULL;
+  GraphNr = -1;
+  ComponentNr = NULL;
+  IsCyclic = false;
+  SeparationVertex = false;
+  LowpointNr = -1;
+  AdaptiveOrder = 0;
+  MaxOrder = false;
   type = pointer->type;  // copy element of atom
   x.CopyVector(&pointer->x); // copy coordination
@@ -54 +60 @@
 atom::~atom()
 {
-  Free((void **)&Name, "atom::~atom: *Name");
   Free((void **)&ComponentNr, "atom::~atom: *ComponentNr");
 };

src/bond.cpp

@@ -5 +5 @@
  */
 
-#include "molecules.hpp"
+#include "bond.hpp"
 
 

src/boundary.cpp

@@ -1 +1 @@
 #include "boundary.hpp"
-#include "linkedcell.hpp"
-#include "molecules.hpp"
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_permutation.h>
-
-#define DEBUG 1
-#define DoSingleStepOutput 0
-#define DoTecplotOutput 1
-#define DoRaster3DOutput 1
-#define DoVRMLOutput 1
-#define TecplotSuffix ".dat"
-#define Raster3DSuffix ".r3d"
-#define VRMLSUffix ".wrl"
-#define HULLEPSILON 1e-7
-
-// ======================================== 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;
-  if (!lines.empty())
-    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
-  node = NULL;
-}
-;
-
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
-{
-  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
-      << 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); //Taken out, to check whether we can avoid unwanted double adding.
-  Point[1]->AddLine(this); //
-  // clear triangles list
-  TrianglesCount = 0;
-  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
-}
-;
-
-BoundaryLineSet::~BoundaryLineSet()
-{
-  int Numbers[2];
-  Numbers[0] = endpoints[1]->Nr;
-  Numbers[1] = endpoints[0]->Nr;
-  for (int i = 0; i < 2; i++) {
-    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-    // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
-    pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
-    for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
-      if ((*Runner).second == this) {
-        endpoints[i]->lines.erase(Runner);
-        break;
-      }
-    if (endpoints[i]->lines.empty()) {
-      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-      if (endpoints[i] != NULL) {
-        delete(endpoints[i]);
-        endpoints[i] = NULL;
-      } else
-        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
-    } else
-      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
-  }
-  if (!triangles.empty())
-    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
-}
-;
-
-void
-BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
-{
-  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
-      << endl;
-  triangles.insert(TrianglePair(triangle->Nr, triangle));
-  TrianglesCount++;
-}
-;
-
-/** Checks whether we have a common endpoint with given \a *line.
- * \param *line other line to test
- * \return true - common endpoint present, false - not connected
- */
-bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
-{
-  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
-    return true;
-  else
-    return false;
-};
-
-/** Checks whether the adjacent triangles of a baseline are convex or not.
- * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
- * If greater/equal M_PI than we are convex.
- * \param *out output stream for debugging
- * \return true - triangles are convex, false - concave or less than two triangles connected
- */
-bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
-{
-  Vector BaseLineNormal;
-  double angle = 0;
-  // get the two triangles
-  if (TrianglesCount != 2) {
-    *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
-    return false;
-  }
-  // have a normal vector on the base line pointing outwards
-  BaseLineNormal.Zero();
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
-    BaseLineNormal.AddVector(&runner->second->NormalVector);
-  BaseLineNormal.Normalize();
-  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
-    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
-
-  if ((angle - M_PI) > -MYEPSILON)
-    return true;
-  else
-    return false;
-}
-
-/** Checks whether point is any of the two endpoints this line contains.
- * \param *point point to test
- * \return true - point is of the line, false - is not
- */
-bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
-{
-  for(int i=0;i<2;i++)
-    if (point == endpoints[i])
-      return true;
-  return false;
-};
-
-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);
-    if (lines[i]->triangles.empty()) {
-      if (lines[i] != NULL) {
-        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-        delete (lines[i]);
-        lines[i] = NULL;
-      } else
-        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
-    } else
-      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
-  }
-}
-;
-
-/** Calculates the normal vector for this triangle.
- * Is made unique by comparison with \a OtherVector to point in the other direction.
- * \param &OtherVector direction vector to make normal vector unique.
- */
-void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
-{
-  // 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 (NormalVector.Projection(&OtherVector) > 0)
-    NormalVector.Scale(-1.);
-};
-
-/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
- * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
- * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
- * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
- * smaller than the first line.
- * \param *out output stream for debugging
- * \param *MolCenter offset vector of line
- * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
- * \param *Intersection intersection on plane on return
- * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
- */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
-{
-  Vector CrossPoint;
-  Vector helper;
-  int i=0;
-
-  if (Intersection->GetIntersectionWithPlane(out, &NormalVector, &endpoints[0]->node->x, MolCenter, x)) {
-    *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
-    return false;
-  }
-
-  // Calculate cross point between one baseline and the line from the third endpoint to intersection
-  do {
-    CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, &endpoints[i%3]->node->x, &endpoints[(i+1)%3]->node->x, &endpoints[(i+2)%3]->node->x, Intersection);
-    helper.CopyVector(&endpoints[(i+1)%3]->node->x);
-    helper.SubtractVector(&endpoints[i%3]->node->x);
-    i++;
-    if (i>3)
-      break;
-  } while (CrossPoint.NormSquared() < MYEPSILON);
-  if (i>3) {
-    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
-    exit(255);
-  }
-  CrossPoint.SubtractVector(&endpoints[i%3]->node->x);
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
-    return true;
-  } else { // outside!
-    Intersection->Zero();
-    return false;
-  }
-};
-
-/** Checks whether lines is any of the three boundary lines this triangle contains.
- * \param *line line to test
- * \return true - line is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
-{
-  for(int i=0;i<3;i++)
-    if (line == lines[i])
-      return true;
-  return false;
-};
-
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point point to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
-{
-  for(int i=0;i<3;i++)
-    if (point == endpoints[i])
-      return true;
-  return false;
-};
-
-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;
-}
-;
-
-
-// ============================ CandidateForTesselation =============================
-
-CandidateForTesselation::CandidateForTesselation(
-  atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
-) {
-  point = candidate;
-  BaseLine = line;
-  OptCenter.CopyVector(&OptCandidateCenter);
-  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
-}
-
-CandidateForTesselation::~CandidateForTesselation() {
-  point = NULL;
-  BaseLine = NULL;
-}
+
+#include<gsl/gsl_poly.h>
 
 // ========================================== 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;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  Vector helper;
-
-  *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 << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << 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.SubtractVector(MolCenter);
-      ProjectedVector.ProjectOntoPlane(&AxisVector);
-
-      // correct for negative side
-      radius = ProjectedVector.NormSquared();
-      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 << endl;
-      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
-      if (!BoundaryTestPair.second) { // 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 << endl;
-        *out << Verbose(2) << "New vector: " << *Walker << endl;
-        double tmp = ProjectedVector.NormSquared();
-        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
-          BoundaryTestPair.first->second.first = tmp;
-          BoundaryTestPair.first->second.second = Walker;
-          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
-        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper.CopyVector(&Walker->x);
-          helper.SubtractVector(MolCenter);
-          tmp = helper.NormSquared();
-          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
-          helper.SubtractVector(MolCenter);
-          if (helper.NormSquared() < tmp) {
-            BoundaryTestPair.first->second.second = Walker;
-            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
-          } else {
-            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
-          }
-        } else {
-          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << 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;
-    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
-    do { // do as long as we still throw one out per round
-      flag = false;
-      Boundaries::iterator left = BoundaryPoints[axis].end();
-      Boundaries::iterator right = BoundaryPoints[axis].end();
-      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-        // set neighbours correctly
-        if (runner == BoundaryPoints[axis].begin()) {
-          left = BoundaryPoints[axis].end();
-        } else {
-          left = runner;
-        }
-        left--;
-        right = runner;
-        right++;
-        if (right == BoundaryPoints[axis].end()) {
-          right = BoundaryPoints[axis].begin();
-        }
-        // check distance
-
-        // 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.SubtractVector(MolCenter);
-          SideA.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideA: ";
-          //          SideA.Output(out);
-          //          *out << endl;
-
-          SideB.CopyVector(&right->second.second->x);
-          SideB.SubtractVector(MolCenter);
-          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.SubtractVector(MolCenter);
-          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)) < -MYEPSILON) {
-            // throw out point
-            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
-            BoundaryPoints[axis].erase(runner);
-            flag = true;
-          }
-        }
-      }
-    } while (flag);
-  }
-  delete(MolCenter);
-  return BoundaryPoints;
-}
-;
 
 /** Determines greatest diameters of a cluster defined by its convex envelope.
@@ -605 +21 @@
   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;
-    }
+  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;
@@ -735 +148 @@
       for (i=0;i<3;i++) { // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
         *vrmlfile << "\t";
       }
@@ -790 +203 @@
       for (i=0;i<3;i++) {  // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
         *rasterfile << "\t";
       }
@@ -821 +234 @@
     *out << Verbose(2) << "The following triangles were created:";
     int Counter = 1;
-    atom *Walker = NULL;
+    TesselPoint *Walker = NULL;
     for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
       Walker = target->second->node;
       LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
     }
     *tecplot << endl;
@@ -837 +250 @@
   }
 }
+
+
+/** 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;
+  Vector *MolCenter = mol->DetermineCenterOfAll(out);
+  Vector helper;
+
+  *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 << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << 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.SubtractVector(MolCenter);
+      ProjectedVector.ProjectOntoPlane(&AxisVector);
+
+      // correct for negative side
+      radius = ProjectedVector.NormSquared();
+      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 << endl;
+      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
+      if (!BoundaryTestPair.second) { // 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 << endl;
+        *out << Verbose(2) << "New vector: " << *Walker << endl;
+        double tmp = ProjectedVector.NormSquared();
+        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
+          BoundaryTestPair.first->second.first = tmp;
+          BoundaryTestPair.first->second.second = Walker;
+          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
+        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
+          helper.CopyVector(&Walker->x);
+          helper.SubtractVector(MolCenter);
+          tmp = helper.NormSquared();
+          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
+          helper.SubtractVector(MolCenter);
+          if (helper.NormSquared() < tmp) {
+            BoundaryTestPair.first->second.second = Walker;
+            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
+          } else {
+            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
+          }
+        } else {
+          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << 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;
+    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
+    do { // do as long as we still throw one out per round
+      flag = false;
+      Boundaries::iterator left = BoundaryPoints[axis].end();
+      Boundaries::iterator right = BoundaryPoints[axis].end();
+      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        // set neighbours correctly
+        if (runner == BoundaryPoints[axis].begin()) {
+          left = BoundaryPoints[axis].end();
+        } else {
+          left = runner;
+        }
+        left--;
+        right = runner;
+        right++;
+        if (right == BoundaryPoints[axis].end()) {
+          right = BoundaryPoints[axis].begin();
+        }
+        // check distance
+
+        // 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.SubtractVector(MolCenter);
+          SideA.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideA: ";
+          //          SideA.Output(out);
+          //          *out << endl;
+
+          SideB.CopyVector(&right->second.second->x);
+          SideB.SubtractVector(MolCenter);
+          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.SubtractVector(MolCenter);
+          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)) < -MYEPSILON) {
+            // throw out point
+            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+            BoundaryPoints[axis].erase(runner);
+            flag = true;
+          }
+        }
+      }
+    } while (flag);
+  }
+  delete(MolCenter);
+  return BoundaryPoints;
+};
 
 /** Tesselates the convex boundary by finding all boundary points.
@@ -959 +546 @@
       << "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++)
+  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.DistanceSquared(
-          &runner->second->endpoints[1]->node->x));
-      b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
-          &runner->second->endpoints[2]->node->x));
-      c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
-          &runner->second->endpoints[1]->node->x));
+      x.CopyVector(runner->second->endpoints[0]->node->node);
+      x.SubtractVector(runner->second->endpoints[1]->node->node);
+      y.CopyVector(runner->second->endpoints[0]->node->node);
+      y.SubtractVector(runner->second->endpoints[2]->node->node);
+      a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
+      b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
+      c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
       G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + 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));
+      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
+      x.Scale(runner->second->endpoints[1]->node->node->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)
@@ -1227 +808 @@
 };
 
-
-// =========================================================== class TESSELATION ===========================================
-
-/** Constructor of class Tesselation.
- */
-Tesselation::Tesselation()
-{
-  PointsOnBoundaryCount = 0;
-  LinesOnBoundaryCount = 0;
-  TrianglesOnBoundaryCount = 0;
-  TriangleFilesWritten = 0;
-}
-;
-
-/** Constructor of class Tesselation.
- * We have to free all points, lines and triangles.
- */
-Tesselation::~Tesselation()
-{
-  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
-    if (runner->second != NULL) {
-      delete (runner->second);
-      runner->second = NULL;
-    } else
-      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
-  }
-}
-;
-
-/** 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, tmp;
-  Vector PlaneVector, TrialVector;
-  PointMap::iterator A, B, C; // three nodes of the first triangle
-  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
-
-  // with A chosen, take each pair B,C and sort
-  if (A != PointsOnBoundary.end())
-    {
-      B = A;
-      B++;
-      for (; B != PointsOnBoundary.end(); B++)
-        {
-          C = B;
-          C++;
-          for (; C != PointsOnBoundary.end(); C++)
-            {
-              tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
-              distance = tmp * tmp;
-              tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
-              distance += tmp * tmp;
-              tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
-              distance += tmp * tmp;
-              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
-                  PointMap::iterator, PointMap::iterator> (B, C)));
-            }
-        }
-    }
-  //    // 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. pick three baselines forming a triangle
-  // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
-  DistanceMultiMap::iterator baseline = DistanceMMap.begin();
-  for (; baseline != DistanceMMap.end(); baseline++)
-    {
-      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
-      // 2. next, we have to check whether all points reside on only one side of the triangle
-      // 3. construct plane vector
-      PlaneVector.MakeNormalVector(&A->second->node->x,
-          &baseline->second.first->second->node->x,
-          &baseline->second.second->second->node->x);
-      *out << Verbose(2) << "Plane vector of candidate triangle is ";
-      PlaneVector.Output(out);
-      *out << endl;
-      // 4. loop over all points
-      double sign = 0.;
-      PointMap::iterator checker = PointsOnBoundary.begin();
-      for (; checker != PointsOnBoundary.end(); checker++)
-        {
-          // (neglecting A,B,C)
-          if ((checker == A) || (checker == baseline->second.first) || (checker
-              == baseline->second.second))
-            continue;
-          // 4a. project onto plane vector
-          TrialVector.CopyVector(&checker->second->node->x);
-          TrialVector.SubtractVector(&A->second->node->x);
-          distance = TrialVector.Projection(&PlaneVector);
-          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
-            continue;
-          *out << Verbose(3) << "Projection of " << checker->second->node->Name
-              << " yields distance of " << distance << "." << endl;
-          tmp = distance / fabs(distance);
-          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
-          if ((sign != 0) && (tmp != sign))
-            {
-              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
-              *out << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leaves "
-                  << checker->second->node->Name << " outside the convex hull."
-                  << endl;
-              break;
-            }
-          else
-            { // note the sign for later
-              *out << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leave "
-                  << checker->second->node->Name << " inside the convex hull."
-                  << endl;
-              sign = tmp;
-            }
-          // 4d. Check whether the point is inside the triangle (check distance to each node
-          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
-          int innerpoint = 0;
-          if ((tmp < A->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x)) && (tmp
-              < A->second->node->x.DistanceSquared(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x);
-          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
-              &A->second->node->x)) && (tmp
-              < baseline->second.first->second->node->x.DistanceSquared(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.DistanceSquared(
-              &baseline->second.second->second->node->x);
-          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x)) && (tmp
-              < baseline->second.second->second->node->x.DistanceSquared(
-                  &A->second->node->x)))
-            innerpoint++;
-          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
-          if (innerpoint == 3)
-            break;
-        }
-      // 5. come this far, all on same side? Then break 1. loop and construct triangle
-      if (checker == PointsOnBoundary.end())
-        {
-          *out << "Looks like we have a candidate!" << endl;
-          break;
-        }
-    }
-  if (baseline != DistanceMMap.end())
-    {
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = A->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = A->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;
-    }
-  else
-    {
-      *out << Verbose(1) << "No starting triangle found." << endl;
-      exit(255);
-    }
-}
-;
-
-/** 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 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)
- *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
- *    -# 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, molecule *mol)
-{
-  bool flag;
-  PointMap::iterator winner;
-  class BoundaryPointSet *peak = NULL;
-  double SmallestAngle, TempAngle;
-  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *MolCenter = NULL;
-  LineMap::iterator LineChecker[2];
-
-  MolCenter = mol->DetermineCenterOfAll(out);
-  // create a first tesselation with the given BoundaryPoints
-  do {
-    flag = false;
-    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
-      if (baseline->second->TrianglesCount == 1) {
-        // 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;
-
-        // get peak point with respect to this base line's only triangle
-        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
-        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
-        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;
-
-        // prepare some auxiliary vectors
-        Vector BaseLineCenter, BaseLine;
-        BaseLineCenter.CopyVector(&baseline->second->endpoints[0]->node->x);
-        BaseLineCenter.AddVector(&baseline->second->endpoints[1]->node->x);
-        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
-        BaseLine.CopyVector(&baseline->second->endpoints[0]->node->x);
-        BaseLine.SubtractVector(&baseline->second->endpoints[1]->node->x);
-
-        // 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 << endl;
-
-        // normal vector of triangle
-        NormalVector.CopyVector(MolCenter);
-        NormalVector.SubtractVector(&CenterVector);
-        BTS->GetNormalVector(NormalVector);
-        NormalVector.CopyVector(&BTS->NormalVector);
-        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
-
-        // vector in propagation direction (out of triangle)
-        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        PropagationVector.MakeNormalVector(&BaseLine, &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 << endl;
-        winner = PointsOnBoundary.end();
-
-        // loop over all points and calculate angle between normal vector of new and present triangle
-        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) << ":" << endl;
-
-            // first check direction, so that triangles don't intersect
-            VirtualNormalVector.CopyVector(&target->second->node->x);
-            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
-            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
-            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
-            *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
-            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
-              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
-              continue;
-            } else
-              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
-
-            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
-            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()) && (LineChecker[0]->second->TrianglesCount == 2))) {
-              *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;
-            }
-            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
-              *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)
-            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
-              *out << Verbose(4) << "Current target is peak!" << endl;
-              continue;
-            }
-
-            // check for linear dependence
-            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.SubtractVector(&target->second->node->x);
-            helper.CopyVector(&baseline->second->endpoints[1]->node->x);
-            helper.SubtractVector(&target->second->node->x);
-            helper.ProjectOntoPlane(&TempVector);
-            if (fabs(helper.NormSquared()) < MYEPSILON) {
-              *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
-              continue;
-            }
-
-            // in case NOT both were found, 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);
-            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.AddVector(&baseline->second->endpoints[1]->node->x);
-            TempVector.AddVector(&target->second->node->x);
-            TempVector.Scale(1./3.);
-            TempVector.SubtractVector(MolCenter);
-            // make it always point outward
-            if (VirtualNormalVector.Projection(&TempVector) < 0)
-              VirtualNormalVector.Scale(-1.);
-            // calculate angle
-            TempAngle = NormalVector.Angle(&VirtualNormalVector);
-            *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
-            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
-              SmallestAngle = TempAngle;
-              winner = target;
-              *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
-            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
-              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
-              helper.CopyVector(&target->second->node->x);
-              helper.SubtractVector(&BaseLineCenter);
-              helper.ProjectOntoPlane(&BaseLine);
-              // ...the one with the smaller angle is the better candidate
-              TempVector.CopyVector(&target->second->node->x);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              TempAngle = TempVector.Angle(&helper);
-              TempVector.CopyVector(&winner->second->node->x);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              if (TempAngle < TempVector.Angle(&helper)) {
-                TempAngle = NormalVector.Angle(&VirtualNormalVector);
-                SmallestAngle = TempAngle;
-                winner = target;
-                *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
-              } else
-                *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
-            } else
-              *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
-          }
-        } // end of loop over all boundary points
-
-        // 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);
-
-  // exit
-  delete(MolCenter);
-};
-
-/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
- * \param *out output stream for debugging
- * \param *mol molecule structure with atoms
- * \return true - all straddling points insert, false - something went wrong
- */
-bool Tesselation::InsertStraddlingPoints(ofstream *out, molecule *mol)
-{
-  Vector Intersection;
-  atom *Walker = mol->start;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  while (Walker->next != mol->end) {  // we only have to go once through all points, as boundary can become only bigger
-    // get the next triangle
-    BTS = FindClosestTriangleToPoint(out, &Walker->x);
-    if (BTS == NULL) {
-      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
-      return false;
-    }
-    // get the intersection point
-    if (BTS->GetIntersectionInsideTriangle(out, MolCenter, &Walker->x, &Intersection)) {
-      // we have the intersection, check whether in- or outside of boundary
-      if ((MolCenter->DistanceSquared(&Walker->x) - MolCenter->DistanceSquared(&Intersection)) < -MYEPSILON) {
-        // inside, next!
-        *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
-      } else {
-        // outside!
-        *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
-        class BoundaryLineSet *OldLines[3], *NewLines[3];
-        class BoundaryPointSet *OldPoints[3], *NewPoint;
-        // store the three old lines and old points
-        for (int i=0;i<3;i++) {
-          OldLines[i] = BTS->lines[i];
-          OldPoints[i] = BTS->endpoints[i];
-        }
-        // add Walker to boundary points
-        AddPoint(Walker);
-        if (BPS[0] == NULL)
-          NewPoint = BPS[0];
-        else
-          continue;
-        // remove triangle
-        TrianglesOnBoundary.erase(BTS->Nr);
-        // create three new boundary lines
-        for (int i=0;i<3;i++) {
-          BPS[0] = NewPoint;
-          BPS[1] = OldPoints[i];
-          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
-          LinesOnBoundaryCount++;
-        }
-        // create three new triangle with new point
-        for (int i=0;i<3;i++) { // find all baselines
-          BLS[0] = OldLines[i];
-          int n = 1;
-          for (int j=0;j<3;j++) {
-            if (NewLines[j]->IsConnectedTo(BLS[0])) {
-              if (n>2) {
-                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
-                return false;
-              } else
-                BLS[n++] = NewLines[j];
-            }
-          }
-          // create the triangle
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-          TrianglesOnBoundaryCount++;
-        }
-      }
-    } else { // something is wrong with FindClosestTriangleToPoint!
-      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
-      return false;
-    }
-  }
-
-  // exit
-  delete(MolCenter);
-  return true;
-};
-
-/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
- * \param *out output stream for debugging
- * \return true - all baselines were corrected, false - there are still concave pieces
- */
-bool Tesselation::CorrectConcaveBaselines(ofstream *out)
-{
-  class BoundaryTriangleSet *triangle[2];
-  class BoundaryLineSet *OldLines[4], *NewLine;
-  class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal;
-  class BoundaryLineSet *Base = NULL;
-  int OldTriangles[2], OldBaseLine;
-  int i;
-  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
-    Base = baseline->second;
-
-    // check convexity
-    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
-      *out << Verbose(3) << Base << " has two convex triangles." << endl;
-    } else { // not convex!
-      // get the two triangles
-      i=0;
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-        triangle[i++] = runner->second;
-      // gather four endpoints and four lines
-      for (int j=0;j<4;j++)
-        OldLines[j] = NULL;
-      for (int j=0;j<2;j++)
-        OldPoints[j] = NULL;
-      i=0;
-      for (int m=0;m<2;m++) { // go over both triangles
-        for (int j=0;j<3;j++) { // all of their endpoints and baselines
-          if (triangle[m]->lines[j] != Base) // pick not the central baseline
-            OldLines[i++] = triangle[m]->lines[j];
-          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
-            OldPoints[m] = triangle[m]->endpoints[j];
-        }
-      }
-      if (i<4) {
-        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-        return false;
-      }
-      for (int j=0;j<4;j++)
-        if (OldLines[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-          return false;
-        }
-      for (int j=0;j<2;j++)
-        if (OldPoints[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-          return false;
-        }
-
-      // remove triangles
-      for (int j=0;j<2;j++) {
-        OldTriangles[j] = triangle[j]->Nr;
-        TrianglesOnBoundary.erase(OldTriangles[j]);
-        triangle[j] = NULL;
-      }
-
-      // remove baseline
-      OldBaseLine = Base->Nr;
-      LinesOnBoundary.erase(OldBaseLine);
-      Base = NULL;
-
-      // construct new baseline (with same number as old one)
-      BPS[0] = OldPoints[0];
-      BPS[1] = OldPoints[1];
-      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
-      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
-
-      // construct new triangles with flipped baseline
-      i=-1;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
-        i=2;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
-        i=3;
-      if (i!=-1) {
-        BLS[0] = OldLines[0];
-        BLS[1] = OldLines[i];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
-
-        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
-        BLS[1] = OldLines[1];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
-      } else {
-        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
-        return false;
-      }
-    }
-  }
-  return true;
-};
-
-
-/** States whether point is in- or outside of a tesselated surface.
- * \param *pointer point to be checked
- * \return true - is inside, false - is outside
- */
-bool Tesselation::IsInside(Vector *pointer)
-{
-
-  // hier kommt dann Saskias Routine hin...
-
-  return true;
-};
-
-/** Finds the closest triangle to a given point.
- * \param *out output stream for debugging
- * \param *x second endpoint of line
- * \return pointer triangle that is closest, NULL if none was found
- */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
-{
-  class BoundaryTriangleSet *triangle = NULL;
-
-  // hier kommt dann Saskias Routine hin...
-
-  return triangle;
-};
-
-/** Adds an atom to the tesselation::PointsOnBoundary list.
- * \param *Walker atom to add
- */
-void
-Tesselation::AddPoint(atom *Walker)
-{
-  PointTestPair InsertUnique;
-  BPS[0] = new class BoundaryPointSet(Walker);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
-  if (InsertUnique.second) // if new point was not present before, increase counter
-    PointsOnBoundaryCount++;
-  else {
-    delete(BPS[0]);
-    BPS[0] = NULL;
-  }
-}
-;
-
-/** Adds point to Tesselation::PointsOnBoundary if not yet present.
- * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
- * @param Candidate point to add
- * @param n index for this point in Tesselation::TPS array
- */
-void
-Tesselation::AddTrianglePoint(atom* Candidate, int n)
-{
-  PointTestPair InsertUnique;
-  TPS[n] = new class BoundaryPointSet(Candidate);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
-  if (InsertUnique.second) { // if new point was not present before, increase counter
-    PointsOnBoundaryCount++;
-  } else {
-    delete TPS[n];
-    cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
-    TPS[n] = (InsertUnique.first)->second;
-  }
-}
-;
-
-/** Function tries to add line from current Points in BPS to BoundaryLineSet.
- * If successful it raises the line count and inserts the new line into the BLS,
- * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
- * @param *a first endpoint
- * @param *b second endpoint
- * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
-  bool insertNewLine = true;
-
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine;
-    pair<LineMap::iterator,LineMap::iterator> FindPair;
-    FindPair = a->lines.equal_range(b->node->nr);
-
-    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-      // If there is a line with less than two attached triangles, we don't need a new line.
-      if (FindLine->second->TrianglesCount < 2) {
-        insertNewLine = false;
-        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
-
-        BPS[0] = FindLine->second->endpoints[0];
-        BPS[1] = FindLine->second->endpoints[1];
-        BLS[n] = FindLine->second;
-
-        break;
-      }
-    }
-  }
-
-  if (insertNewLine) {
-    AlwaysAddTriangleLine(a, b, n);
-  }
-}
-;
-
-/**
- * Adds lines from each of the current points in the BPS to BoundaryLineSet.
- * Raises the line count and inserts the new line into the BLS.
- *
- * @param *a first endpoint
- * @param *b second endpoint
- * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
-{
-  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
-  BPS[0] = a;
-  BPS[1] = b;
-  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
-  // add line to global map
-  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
-  // increase counter
-  LinesOnBoundaryCount++;
-};
-
-/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
- * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
- */
-void
-Tesselation::AddTriangle()
-{
-  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
-
-  // add triangle to global map
-  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-  TrianglesOnBoundaryCount++;
-
-  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-}
-;
-
-
-double det_get(gsl_matrix *A, int inPlace) {
-  /*
-  inPlace = 1 => A is replaced with the LU decomposed copy.
-  inPlace = 0 => A is retained, and a copy is used for LU.
-  */
-
-  double det;
-  int signum;
-  gsl_permutation *p = gsl_permutation_alloc(A->size1);
-  gsl_matrix *tmpA;
-
-  if (inPlace)
-  tmpA = A;
-  else {
-  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
-  gsl_matrix_memcpy(tmpA , A);
-  }
-
-
-  gsl_linalg_LU_decomp(tmpA , p , &signum);
-  det = gsl_linalg_LU_det(tmpA , signum);
-  gsl_permutation_free(p);
-  if (! inPlace)
-  gsl_matrix_free(tmpA);
-
-  return det;
-};
-
-void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
-{
-  gsl_matrix *A = gsl_matrix_calloc(3,3);
-  double m11, m12, m13, m14;
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
-  }
-  m11 = det_get(A, 1);
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
-  }
-  m12 = det_get(A, 1);
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
-  }
-  m13 = det_get(A, 1);
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, b.x[i]);
-  }
-  m14 = det_get(A, 1);
-
-  if (fabs(m11) < MYEPSILON)
-    cerr << "ERROR: three points are colinear." << endl;
-
-  center->x[0] =  0.5 * m12/ m11;
-  center->x[1] = -0.5 * m13/ m11;
-  center->x[2] =  0.5 * m14/ m11;
-
-  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
-    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
-
-  gsl_matrix_free(A);
-};
-
-
-
-/**
- * Function returns center of sphere with RADIUS, which rests on points a, b, c
- * @param Center this vector will be used for return
- * @param a vector first point of triangle
- * @param b vector second point of triangle
- * @param c vector third point of triangle
- * @param *Umkreismittelpunkt new cneter point of circumference
- * @param Direction vector indicates up/down
- * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
- * @param Halfplaneindicator double indicates whether Direction is up or down
- * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
- * @param alpha double angle at a
- * @param beta double, angle at b
- * @param gamma, double, angle at c
- * @param Radius, double
- * @param Umkreisradius double radius of circumscribing circle
- */
-void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
-    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
-{
-  Vector TempNormal, helper;
-  double Restradius;
-  Vector OtherCenter;
-  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
-  Center->Zero();
-  helper.CopyVector(&a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(&b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(&c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
-  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
-  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-  NewUmkreismittelpunkt->CopyVector(Center);
-  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
-  // Here we calculated center of circumscribing circle, using barycentric coordinates
-  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
-
-  TempNormal.CopyVector(&a);
-  TempNormal.SubtractVector(&b);
-  helper.CopyVector(&a);
-  helper.SubtractVector(&c);
-  TempNormal.VectorProduct(&helper);
-  if (fabs(HalfplaneIndicator) < MYEPSILON)
-    {
-      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
-        {
-          TempNormal.Scale(-1);
-        }
-    }
-  else
-    {
-      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
-        {
-          TempNormal.Scale(-1);
-        }
-    }
-
-  TempNormal.Normalize();
-  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
-  TempNormal.Scale(Restradius);
-  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
-
-  Center->AddVector(&TempNormal);
-  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
-  get_sphere(&OtherCenter, a, b, c, RADIUS);
-  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
-  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
-};
-
-
-/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
- * \param *Center new center on return
- * \param *a first point
- * \param *b second point
- * \param *c third point
- */
-void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
-{
-  Vector helper;
-  double alpha, beta, gamma;
-  Vector SideA, SideB, SideC;
-  SideA.CopyVector(b);
-  SideA.SubtractVector(c);
-  SideB.CopyVector(c);
-  SideB.SubtractVector(a);
-  SideC.CopyVector(a);
-  SideC.SubtractVector(b);
-  alpha = M_PI - SideB.Angle(&SideC);
-  beta = M_PI - SideC.Angle(&SideA);
-  gamma = M_PI - SideA.Angle(&SideB);
-  //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
-  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
-    cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
-
-  Center->Zero();
-  helper.CopyVector(a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
-  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-};
-
-/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
- * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
- * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
- * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
- * \param CircleCenter Center of the parameter circle
- * \param CirclePlaneNormal normal vector to plane of the parameter circle
- * \param CircleRadius radius of the parameter circle
- * \param NewSphereCenter new center of a circumcircle
- * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
- * \param NormalVector normal vector
- * \param SearchDirection search direction to make angle unique on return.
- * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
- */
-double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
-{
-  Vector helper;
-  double radius, alpha;
-
-  helper.CopyVector(&NewSphereCenter);
-  // test whether new center is on the parameter circle's plane
-  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
-    helper.ProjectOntoPlane(&CirclePlaneNormal);
-  }
-  radius = helper.ScalarProduct(&helper);
-  // test whether the new center vector has length of CircleRadius
-  if (fabs(radius - CircleRadius) > HULLEPSILON)
-    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-  alpha = helper.Angle(&OldSphereCenter);
-  // make the angle unique by checking the halfplanes/search direction
-  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
-    alpha = 2.*M_PI - alpha;
-  //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&OldSphereCenter);
-  helper.ProjectOntoPlane(&NormalVector);
-  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
-  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
-    //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
-    return alpha;
-  } else {
-    //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
-    return 2.*M_PI;
-  }
-};
-
-
-/** Checks whether the triangle consisting of the three atoms is already present.
- * Searches for the points in Tesselation::PointsOnBoundary and checks their
- * lines. If any of the three edges already has two triangles attached, false is
- * returned.
- * \param *out output stream for debugging
- * \param *Candidates endpoints of the triangle candidate
- * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
- *                 triangles exist which is the maximum for three points
- */
-int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
-  LineMap::iterator FindLine;
-  PointMap::iterator FindPoint;
-  TriangleMap::iterator FindTriangle;
-  int adjacentTriangleCount = 0;
-  class BoundaryPointSet *Points[3];
-
-  //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
-  // builds a triangle point set (Points) of the end points
-  for (int i = 0; i < 3; i++) {
-    FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      Points[i] = FindPoint->second;
-    } else {
-      Points[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (Points[i] != NULL) {
-      for (int j = i; j < 3; j++) {
-        if (Points[j] != NULL) {
-          FindLine = Points[i]->lines.find(Points[j]->node->nr);
-          if (FindLine != Points[i]->lines.end()) {
-            for (; FindLine->first == Points[j]->node->nr; FindLine++) {
-              FindTriangle = FindLine->second->triangles.begin();
-              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == Points[0])
-                    || (FindTriangle->second->endpoints[0] == Points[1])
-                    || (FindTriangle->second->endpoints[0] == Points[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == Points[0])
-                    || (FindTriangle->second->endpoints[1] == Points[1])
-                    || (FindTriangle->second->endpoints[1] == Points[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == Points[0])
-                    || (FindTriangle->second->endpoints[2] == Points[1])
-                    || (FindTriangle->second->endpoints[2] == Points[2])
-                  )
-                ) {
-                  adjacentTriangleCount++;
-                }
-              }
-            }
-            // Only one of the triangle lines must be considered for the triangle count.
-            *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-            return adjacentTriangleCount;
-
-          }
-        }
-      }
-    }
-  }
-
-  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-  return adjacentTriangleCount;
-};
-
-/** This recursive function finds a third point, to form a triangle with two given ones.
- * Note that this function is for the starting triangle.
- * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
- * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
- * the center of the sphere is still fixed up to a single parameter. The band of possible values
- * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
- * us the "null" on this circle, the new center of the candidate point will be some way along this
- * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
- * by the normal vector of the base triangle that always points outwards by construction.
- * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
- * We construct the normal vector that defines the plane this circle lies in, it is just in the
- * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
- * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
- * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
- * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
- * both.
- * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
- * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
- * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
- * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
- * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
- * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
- * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
- * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
- * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
- * @param BaseLine BoundaryLineSet with the current base line
- * @param ThirdNode third atom to avoid in search
- * @param candidates list of equally good candidates to return
- * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
- * @param RADIUS radius of sphere
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-void Find_third_point_for_Tesselation(
-  Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter,
-  class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates,
-  double *ShortestAngle, const double RADIUS, LinkedCell *LC
-) {
-  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
-  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
-  Vector SphereCenter;
-  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
-  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
-  Vector NewNormalVector;   // normal vector of the Candidate's triangle
-  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
-  LinkedAtoms *List = NULL;
-  double CircleRadius; // radius of this circle
-  double radius;
-  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-  atom *Candidate = NULL;
-  CandidateForTesselation *optCandidate = NULL;
-
-  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
-
-  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
-
-  // construct center of circle
-  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
-  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
-  CircleCenter.Scale(0.5);
-
-  // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
-  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
-
-  // calculate squared radius atom *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
-    CirclePlaneNormal.Normalize();
-    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-
-    // test whether old center is on the band's plane
-    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
-    }
-    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
-    if (fabs(radius - CircleRadius) < HULLEPSILON) {
-
-      // check SearchDirection
-      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
-        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
-      }
-
-      // get cell for the starting atom
-      if (LC->SetIndexToVector(&CircleCenter)) {
-        for(int i=0;i<NDIM;i++) // store indices of this cell
-        N[i] = LC->n[i];
-        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-      } else {
-        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
-        return;
-      }
-      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-      //cout << Verbose(2) << "LC Intervals:";
-      for (int i=0;i<NDIM;i++) {
-        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-      }
-      //cout << endl;
-      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-            List = LC->GetCurrentCell();
-            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-            if (List != NULL) {
-              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-                Candidate = (*Runner);
-
-                // check for three unique points
-                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
-                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
-
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-
-                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
-                  ) {
-                    helper.CopyVector(&NewNormalVector);
-                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
-                    if (radius < RADIUS*RADIUS) {
-                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
-                      NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
-                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-
-                      // OtherNewSphereCenter is created by the same vector just in the other direction
-                      helper.Scale(-1.);
-                      OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
-                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
-
-                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                      alpha = min(alpha, Otheralpha);
-                      // if there is a better candidate, drop the current list and add the new candidate
-                      // otherwise ignore the new candidate and keep the list
-                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
-                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
-                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
-                        } else {
-                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
-                        }
-                        // if there is an equal candidate, add it to the list without clearing the list
-                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
-                          candidates->push_back(optCandidate);
-                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
-                        } else {
-                          // remove all candidates from the list and then the list itself
-                          class CandidateForTesselation *remover = NULL;
-                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
-                            remover = *it;
-                            delete(remover);
-                          }
-                          candidates->clear();
-                          candidates->push_back(optCandidate);
-                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
-                        }
-                        *ShortestAngle = alpha;
-                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
-                      } else {
-                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
-                          //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
-                        } else {
-                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
-                        }
-                      }
-
-                    } else {
-                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
-                    }
-                  } else {
-                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
-                  }
-                } else {
-                  if (ThirdNode != NULL) {
-                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
-                  } else {
-                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
-                  }
-                }
-              }
-            }
-          }
-    } else {
-      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-    }
-  } else {
-    if (ThirdNode != NULL)
-      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
-    else
-      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
-  }
-
-  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
-  if (candidates->size() > 1) {
-    candidates->unique();
-    candidates->sort(sortCandidates);
-  }
-
-  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
-};
-
-struct Intersection {
-  Vector x1;
-  Vector x2;
-  Vector x3;
-  Vector x4;
-};
-
-/**
- * Intersection calculation function.
- *
- * @param x to find the result for
- * @param function parameter
- */
-double MinIntersectDistance(const gsl_vector * x, void *params) {
-  double retval = 0;
-  struct Intersection *I = (struct Intersection *)params;
-  Vector intersection;
-  Vector SideA,SideB,HeightA, HeightB;
-  for (int i=0;i<NDIM;i++)
-    intersection.x[i] = gsl_vector_get(x, i);
-
-  SideA.CopyVector(&(I->x1));
-  SideA.SubtractVector(&I->x2);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&I->x1);
-  HeightA.ProjectOntoPlane(&SideA);
-
-  SideB.CopyVector(&I->x3);
-  SideB.SubtractVector(&I->x4);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&I->x3);
-  HeightB.ProjectOntoPlane(&SideB);
-
-  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
-  //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
-
-  return retval;
-};
-
-
-/**
- * Calculates whether there is an intersection between two lines. The first line
- * always goes through point 1 and point 2 and the second line is given by the
- * connection between point 4 and point 5.
- *
- * @param point 1 of line 1
- * @param point 2 of line 1
- * @param point 1 of line 2
- * @param point 2 of line 2
- *
- * @return true if there is an intersection between the given lines, false otherwise
- */
-bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
-  bool result;
-
-  struct Intersection par;
-    par.x1.CopyVector(&point1);
-    par.x2.CopyVector(&point2);
-    par.x3.CopyVector(&point3);
-    par.x4.CopyVector(&point4);
-
-    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
-    gsl_multimin_fminimizer *s = NULL;
-    gsl_vector *ss, *x;
-    gsl_multimin_function minex_func;
-
-    size_t iter = 0;
-    int status;
-    double size;
-
-    /* Starting point */
-    x = gsl_vector_alloc(NDIM);
-    gsl_vector_set(x, 0, point1.x[0]);
-    gsl_vector_set(x, 1, point1.x[1]);
-    gsl_vector_set(x, 2, point1.x[2]);
-
-    /* Set initial step sizes to 1 */
-    ss = gsl_vector_alloc(NDIM);
-    gsl_vector_set_all(ss, 1.0);
-
-    /* Initialize method and iterate */
-    minex_func.n = NDIM;
-    minex_func.f = &MinIntersectDistance;
-    minex_func.params = (void *)&par;
-
-    s = gsl_multimin_fminimizer_alloc(T, NDIM);
-    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
-
-    do {
-        iter++;
-        status = gsl_multimin_fminimizer_iterate(s);
-
-        if (status) {
-          break;
-        }
-
-        size = gsl_multimin_fminimizer_size(s);
-        status = gsl_multimin_test_size(size, 1e-2);
-
-        if (status == GSL_SUCCESS) {
-          cout << Verbose(2) << "converged to minimum" <<  endl;
-        }
-    } while (status == GSL_CONTINUE && iter < 100);
-
-    // check whether intersection is in between or not
-  Vector intersection, SideA, SideB, HeightA, HeightB;
-  double t1, t2;
-  for (int i = 0; i < NDIM; i++) {
-    intersection.x[i] = gsl_vector_get(s->x, i);
-  }
-
-  SideA.CopyVector(&par.x2);
-  SideA.SubtractVector(&par.x1);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&par.x1);
-
-  t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
-
-  SideB.CopyVector(&par.x4);
-  SideB.SubtractVector(&par.x3);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&par.x3);
-
-  t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
-
-  cout << Verbose(2) << "Intersection " << intersection << " is at "
-    << t1 << " for (" << point1 << "," << point2 << ") and at "
-    << t2 << " for (" << point3 << "," << point4 << "): ";
-
-  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
-    cout << "true intersection." << endl;
-    result = true;
-  } else {
-    cout << "intersection out of region of interest." << endl;
-    result = false;
-  }
-
-  // free minimizer stuff
-    gsl_vector_free(x);
-    gsl_vector_free(ss);
-    gsl_multimin_fminimizer_free(s);
-
-  return result;
-}
-
-/** Finds the second point of starting triangle.
- * \param *a first atom
- * \param *Candidate pointer to candidate atom on return
- * \param Oben vector indicating the outside
- * \param Opt_Candidate reference to recommended candidate on return
- * \param Storage[3] array storing angles and other candidate information
- * \param RADIUS radius of virtual sphere
- * \param *LC LinkedCell structure with neighbouring atoms
- */
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
-{
-  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
-  Vector AngleCheck;
-  double norm = -1., angle;
-  LinkedAtoms *List = NULL;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
-    for(int i=0;i<NDIM;i++) // store indices of this cell
-      N[i] = LC->n[i];
-  } else {
-    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
-    return;
-  }
-  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-  cout << Verbose(3) << "LC Intervals from [";
-  for (int i=0;i<NDIM;i++) {
-  cout << " " << N[i] << "<->" << LC->N[i];
-  }
-  cout << "] :";
-  for (int i=0;i<NDIM;i++) {
-    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-  }
-  cout << endl;
-
-
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            Candidate = (*Runner);
-            // check if we only have one unique point yet ...
-            if (a != Candidate) {
-              // Calculate center of the circle with radius RADIUS through points a and Candidate
-              Vector OrthogonalizedOben, a_Candidate, Center;
-              double distance, scaleFactor;
-
-              OrthogonalizedOben.CopyVector(&Oben);
-              a_Candidate.CopyVector(&(a->x));
-              a_Candidate.SubtractVector(&(Candidate->x));
-              OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
-              OrthogonalizedOben.Normalize();
-              distance = 0.5 * a_Candidate.Norm();
-              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
-              OrthogonalizedOben.Scale(scaleFactor);
-
-              Center.CopyVector(&(Candidate->x));
-              Center.AddVector(&(a->x));
-              Center.Scale(0.5);
-              Center.AddVector(&OrthogonalizedOben);
-
-              AngleCheck.CopyVector(&Center);
-              AngleCheck.SubtractVector(&(a->x));
-              norm = a_Candidate.Norm();
-              // second point shall have smallest angle with respect to Oben vector
-              if (norm < RADIUS*2.) {
-                angle = AngleCheck.Angle(&Oben);
-                if (angle < Storage[0]) {
-                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
-                  Opt_Candidate = Candidate;
-                  Storage[0] = angle;
-                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-                } else {
-                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
-                }
-              } else {
-                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
-              }
-            } else {
-              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
-            }
-          }
-        } else {
-          cout << Verbose(3) << "Linked cell list is empty." << endl;
-        }
-      }
-  cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
-};
-
-/** Finds the starting triangle for find_non_convex_border().
- * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
- * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
- * point are called.
- * \param RADIUS radius of virtual rolling sphere
- * \param *LC LinkedCell structure with neighbouring atoms
- */
-void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
-{
-  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
-  int i = 0;
-  LinkedAtoms *List = NULL;
-  atom* FirstPoint = NULL;
-  atom* SecondPoint = NULL;
-  atom* MaxAtom[NDIM];
-  double max_coordinate[NDIM];
-  Vector Oben;
-  Vector helper;
-  Vector Chord;
-  Vector SearchDirection;
-
-  Oben.Zero();
-
-  for (i = 0; i < 3; i++) {
-    MaxAtom[i] = NULL;
-    max_coordinate[i] = -1;
-  }
-
-  // 1. searching topmost atom with respect to each axis
-  for (int i=0;i<NDIM;i++) { // each axis
-    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
-    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
-      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
-        List = LC->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if ((*Runner)->x.x[i] > max_coordinate[i]) {
-              cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl;
-              max_coordinate[i] = (*Runner)->x.x[i];
-              MaxAtom[i] = (*Runner);
-            }
-          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-  }
-
-  cout << Verbose(2) << "Found maximum coordinates: ";
-  for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *MaxAtom[i] << "\t";
-  cout << endl;
-
-  BTS = NULL;
-  CandidateList *Opt_Candidates = new CandidateList();
-  for (int k=0;k<NDIM;k++) {
-    Oben.x[k] = 1.;
-    FirstPoint = MaxAtom[k];
-    cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
-
-    double ShortestAngle;
-    atom* Opt_Candidate = NULL;
-    ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-
-    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
-    SecondPoint = Opt_Candidate;
-    if (SecondPoint == NULL)  // have we found a second point?
-      continue;
-    else
-      cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-
-    helper.CopyVector(&(FirstPoint->x));
-    helper.SubtractVector(&(SecondPoint->x));
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
-    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-
-    Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-    Chord.SubtractVector(&(SecondPoint->x));
-    double radius = Chord.ScalarProduct(&Chord);
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
-    helper.CopyVector(&Oben);
-    helper.Scale(CircleRadius);
-    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-
-    // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
-
-    // adding point 1 and point 2 and the line between them
-    AddTrianglePoint(FirstPoint, 0);
-    AddTrianglePoint(SecondPoint, 1);
-    AddTriangleLine(TPS[0], TPS[1], 0);
-
-    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    Find_third_point_for_Tesselation(
-      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
-    );
-    cout << Verbose(1) << "List of third Points is ";
-    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-        cout << " " << *(*it)->point;
-    }
-    cout << endl;
-
-    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-      // add third triangle point
-      AddTrianglePoint((*it)->point, 2);
-      // add the second and third line
-      AddTriangleLine(TPS[1], TPS[2], 1);
-      AddTriangleLine(TPS[0], TPS[2], 2);
-      // ... and triangles to the Maps of the Tesselation class
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
-      // ... and calculate its normal vector (with correct orientation)
-      (*it)->OptCenter.Scale(-1.);
-      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
-      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
-      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
-      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
-
-      // if we do not reach the end with the next step of iteration, we need to setup a new first line
-      if (it != Opt_Candidates->end()--) {
-        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
-        SecondPoint = (*it)->point;
-        // adding point 1 and point 2 and the line between them
-        AddTrianglePoint(FirstPoint, 0);
-        AddTrianglePoint(SecondPoint, 1);
-        AddTriangleLine(TPS[0], TPS[1], 0);
-      }
-      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
-    }
-    if (BTS != NULL) // we have created one starting triangle
-      break;
-    else {
-      // remove all candidates from the list and then the list itself
-      class CandidateForTesselation *remover = NULL;
-      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-        remover = *it;
-        delete(remover);
-      }
-      Opt_Candidates->clear();
-    }
-  }
-
-  // remove all candidates from the list and then the list itself
-  class CandidateForTesselation *remover = NULL;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    remover = *it;
-    delete(remover);
-  }
-  delete(Opt_Candidates);
-  cout << Verbose(1) << "End of Find_starting_triangle\n";
-};
-
-/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
- * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
- * make it bigger (i.e. closing one (the baseline) and opening two new ones).
- * \param TPS[3] nodes of the triangle
- * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
- */
-bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
-{
-  bool result = false;
-  int counter = 0;
-
-  // check all three points
-  for (int i=0;i<3;i++)
-    for (int j=i+1; j<3; j++) {
-      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
-        LineMap::iterator FindLine;
-        pair<LineMap::iterator,LineMap::iterator> FindPair;
-        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
-        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-          // If there is a line with less than two attached triangles, we don't need a new line.
-          if (FindLine->second->TrianglesCount < 2) {
-            counter++;
-            break;  // increase counter only once per edge
-          }
-        }
-      } else { // no line
-        cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
-        result = true;
-      }
-    }
-  if (counter > 1) {
-    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
-    result = true;
-  }
-  return result;
-};
-
-
-/** This function finds a triangle to a line, adjacent to an existing one.
- * @param out output stream for debugging
- * @param *mol molecule with Atom's and Bond's
- * @param Line current baseline to search from
- * @param T current triangle which \a Line is edge of
- * @param RADIUS radius of the rolling ball
- * @param N number of found triangles
- * @param *filename filename base for intermediate envelopes
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-bool Tesselation::Find_next_suitable_triangle(ofstream *out,
-    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
-{
-  cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
-  ofstream *tempstream = NULL;
-  char NumberName[255];
-  bool result = true;
-  CandidateList *Opt_Candidates = new CandidateList();
-
-  Vector CircleCenter;
-  Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
-  Vector SearchDirection;
-  Vector helper;
-  atom *ThirdNode = NULL;
-  LineMap::iterator testline;
-  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-  double radius, CircleRadius;
-
-  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-  for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
-      ThirdNode = T.endpoints[i]->node;
-
-  // construct center of circle
-  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
-  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
-  CircleCenter.Scale(0.5);
-
-  // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
-  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
-
-  // calculate squared radius of circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
-    CirclePlaneNormal.Normalize();
-    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-
-    // construct old center
-    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
-    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(&Line.endpoints[0]->node->x);
-    helper.SubtractVector(&ThirdNode->x);
-    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
-      SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
-    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
-      // rotated the wrong way!
-      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
-    }
-
-    // add third point
-    Find_third_point_for_Tesselation(
-      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
-      &ShortestAngle, RADIUS, LC
-    );
-
-  } else {
-    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
-  }
-
-  if (Opt_Candidates->begin() == Opt_Candidates->end()) {
-    cerr << "WARNING: Could not find a suitable candidate." << endl;
-    return false;
-  }
-  cout << Verbose(1) << "Third Points are ";
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    cout << " " << *(*it)->point;
-  }
-  cout << endl;
-
-  BoundaryLineSet *BaseRay = &Line;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    cout << Verbose(1) << " Third point candidate is " << *(*it)->point
-    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
-    cout << Verbose(1) << " Baseline is " << *BaseRay << endl;
-
-    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-    atom *AtomCandidates[3];
-    AtomCandidates[0] = (*it)->point;
-    AtomCandidates[1] = BaseRay->endpoints[0]->node;
-    AtomCandidates[2] = BaseRay->endpoints[1]->node;
-    int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
-
-    BTS = NULL;
-    // If there is no triangle, add it regularly.
-    if (existentTrianglesCount == 0) {
-      AddTrianglePoint((*it)->point, 0);
-      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
-
-      AddTriangleLine(TPS[0], TPS[1], 0);
-      AddTriangleLine(TPS[0], TPS[2], 1);
-      AddTriangleLine(TPS[1], TPS[2], 2);
-
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
-      (*it)->OptCenter.Scale(-1.);
-      BTS->GetNormalVector((*it)->OptCenter);
-      (*it)->OptCenter.Scale(-1.);
-
-      cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-        << " for this triangle ... " << endl;
-      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
-    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
-        AddTrianglePoint((*it)->point, 0);
-        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
-
-        // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
-        // i.e. at least one of the three lines must be present with TriangleCount <= 1
-        if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
-          AddTriangleLine(TPS[0], TPS[1], 0);
-          AddTriangleLine(TPS[0], TPS[2], 1);
-          AddTriangleLine(TPS[1], TPS[2], 2);
-
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTriangle();
-
-          (*it)->OtherOptCenter.Scale(-1.);
-          BTS->GetNormalVector((*it)->OtherOptCenter);
-          (*it)->OtherOptCenter.Scale(-1.);
-
-          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-          << " for this triangle ... " << endl;
-          cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
-        } else {
-          cout << Verbose(1) << "WARNING: This triangle consisting of ";
-          cout << *(*it)->point << ", ";
-          cout << *BaseRay->endpoints[0]->node << " and ";
-          cout << *BaseRay->endpoints[1]->node << " ";
-          cout << "exists and is not added, as it does not seem helpful!" << endl;
-          result = false;
-        }
-    } else {
-      cout << Verbose(1) << "This triangle consisting of ";
-      cout << *(*it)->point << ", ";
-      cout << *BaseRay->endpoints[0]->node << " and ";
-      cout << *BaseRay->endpoints[1]->node << " ";
-      cout << "is invalid!" << endl;
-      result = false;
-    }
-
-    if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
-      if (DoTecplotOutput) {
-        string NameofTempFile(tempbasename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(TecplotSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-
-      if (DoRaster3DOutput) {
-        string NameofTempFile(tempbasename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(Raster3DSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, this, mol);
-        // include the current position of the virtual sphere in the temporary raster3d file
-        // make the circumsphere's center absolute again
-        helper.CopyVector(&BaseRay->endpoints[0]->node->x);
-        helper.AddVector(&BaseRay->endpoints[1]->node->x);
-        helper.Scale(0.5);
-        (*it)->OptCenter.AddVector(&helper);
-        Vector *center = mol->DetermineCenterOfAll(out);
-        (*it)->OptCenter.SubtractVector(center);
-        delete(center);
-        // and add to file plus translucency object
-        *tempstream << "# current virtual sphere\n";
-        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-          << "\t" << RADIUS << "\t1 0 0\n";
-        *tempstream << "9\n  terminating special property\n";
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-      if (DoTecplotOutput || DoRaster3DOutput)
-        TriangleFilesWritten++;
-    }
-
-    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
-    BaseRay = BLS[0];
-  }
-
-  // remove all candidates from the list and then the list itself
-  class CandidateForTesselation *remover = NULL;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    remover = *it;
-    delete(remover);
-  }
-  delete(Opt_Candidates);
-  cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
-  return result;
-};
-
-/**
- * Sort function for the candidate list.
- */
-bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
-  Vector BaseLineVector, OrthogonalVector, helper;
-  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
-    //return false;
-    exit(1);
-  }
-  // create baseline vector
-  BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
-  BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  BaseLineVector.Normalize();
-
-  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
-  helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  helper.SubtractVector(&(candidate1->point->x));
-  OrthogonalVector.CopyVector(&helper);
-  helper.VectorProduct(&BaseLineVector);
-  OrthogonalVector.SubtractVector(&helper);
-  OrthogonalVector.Normalize();
-
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(&(candidate1->point->x));
-  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  double phi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-  helper.CopyVector(&(candidate2->point->x));
-  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  double psi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    psi = 2.*M_PI - psi;
-  }
-
-  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
-  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
-
-  // return comparison
-  return phi < psi;
-}
-
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
  * \param *out output stream for debugging
@@ -3176 +821 @@
   bool freeTess = false;
   bool freeLC = false;
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  int TriangleFilesWritten = 0;
+
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
   if (Tess == NULL) {
@@ -3193 +842 @@
   }
 
-  Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
+  Tess->Find_starting_triangle(out, RADIUS, LCList);
 
   baseline = Tess->LinesOnBoundary.begin();
   while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
     if (baseline->second->TrianglesCount == 1) {
-      failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
+      failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
@@ -3214 +863 @@
       flag = false;
     }
-  }
+
+    // write temporary envelope
+    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
+      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
+      if (DoTecplotOutput) {
+        string NameofTempFile(filename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(TecplotSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+
+      if (DoRaster3DOutput) {
+        string NameofTempFile(filename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(Raster3DSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_raster3d_file(out, tempstream, Tess, mol);
+//        // include the current position of the virtual sphere in the temporary raster3d file
+//        // make the circumsphere's center absolute again
+//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
+//        helper.AddVector(BaseRay->endpoints[1]->node->node);
+//        helper.Scale(0.5);
+//        (*it)->OptCenter.AddVector(&helper);
+//        Vector *center = mol->DetermineCenterOfAll(out);
+//        (*it)->OptCenter.SubtractVector(center);
+//        delete(center);
+//        // and add to file plus translucency object
+//        *tempstream << "# current virtual sphere\n";
+//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
+//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
+//          << "\t" << RADIUS << "\t1 0 0\n";
+//        *tempstream << "9\n  terminating special property\n";
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+      if (DoTecplotOutput || DoRaster3DOutput)
+        TriangleFilesWritten++;
+    }
+  }
+
+  // write final envelope
   if (filename != 0) {
     *out << Verbose(1) << "Writing final tecplot file\n";

src/boundary.hpp

@@ -1 +1 @@
 #ifndef BOUNDARY_HPP_
 #define BOUNDARY_HPP_
-
-class BoundaryPointSet;
-class BoundaryLineSet;
-class BoundaryTriangleSet;
-class CandidateForTesselation;
 
 // include config.h
@@ -17 +12 @@
 #include <deque>
 
-#include <gsl/gsl_poly.h>
-
 #include "linkedcell.hpp"
 #include "molecules.hpp"
+#include "tesselation.hpp"
+#include "tesselationhelpers.hpp"
 
-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);
-    bool IsConnectedTo(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-    bool CheckConvexityCriterion(ofstream *out);
-
-    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);
-    bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
-    bool ContainsBoundaryLine(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-
-    class BoundaryPointSet *endpoints[3];
-    class BoundaryLineSet *lines[3];
-    Vector NormalVector;
-    int Nr;
-};
-
-
-class CandidateForTesselation {
-  public :
-  CandidateForTesselation(atom* candidate, BoundaryLineSet* currentBaseLine, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
-  ~CandidateForTesselation();
-  atom *point;
-  BoundaryLineSet *BaseLine;
-  Vector OptCenter;
-  Vector OtherOptCenter;
-};
-
-
-class Tesselation {
-  public:
-
-    Tesselation();
-    ~Tesselation();
-
-    void TesselateOnBoundary(ofstream *out, molecule *mol);
-    void GuessStartingTriangle(ofstream *out);
-    void AddPoint(atom * Walker);
-    void AddTrianglePoint(atom* Candidate, int n);
-    void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AddTriangle();
-    void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC);
-    bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC);
-    int CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]);
-    void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol);
-    class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, Vector *x);
-    bool IsInside(Vector *pointer);
-    bool InsertStraddlingPoints(ofstream *out, molecule *mol);
-    bool CorrectConcaveBaselines(ofstream *out);
-
-    PointMap PointsOnBoundary;
-    LineMap LinesOnBoundary;
-    TriangleMap TrianglesOnBoundary;
-    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
-    class BoundaryPointSet *BPS[2];
-    class BoundaryLineSet *BLS[3];
-    class BoundaryTriangleSet *BTS;
-    int PointsOnBoundaryCount;
-    int LinesOnBoundaryCount;
-    int TrianglesOnBoundaryCount;
-    int TriangleFilesWritten;
-};
-
-
-ostream & operator << (ostream &ost, BoundaryPointSet &a);
-ostream & operator << (ostream &ost, BoundaryLineSet &a);
-ostream & operator << (ostream &ost, BoundaryTriangleSet &a);
-
+#define DEBUG 1
+#define DoSingleStepOutput 0
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 1
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
 
 double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration);
@@ -141 +33 @@
 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
-bool Choose_preferable_third_point(atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, atom *ThirdNode, Tesselation *Tess);
-bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4);
-bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2);
+Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
 
 #endif /*BOUNDARY_HPP_*/

src/builder.cpp

@@ -1229 +1229 @@
   }
   molecules->SimpleMultiAdd(mol, src, N);
-  delete[](src);
+  delete(src);
+
   // ... and translate back
   for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {

src/config.cpp

@@ -276 +276 @@
 
 /** Readin of Thermostat related values from parameter file.
- * \param *source parameter file
- */
-void config::InitThermostats(ifstream *source)
+ * \param *fb file buffer containing the config file
+ */
+void config::InitThermostats(class ConfigFileBuffer *fb)
 {
   char *thermo = MallocString(12, "IonsInitRead: thermo");
@@ -284 +284 @@
 
   // read desired Thermostat from file along with needed additional parameters
-  if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
+  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
     if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
       if (ThermostatImplemented[0] == 1) {
@@ -295 +295 @@
       if (ThermostatImplemented[1] == 1) {
         Thermostat = Woodcock;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
       } else {
         cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
@@ -303 +303 @@
       if (ThermostatImplemented[2] == 1) {
         Thermostat = Gaussian;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
       } else {
         cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
@@ -311 +311 @@
       if (ThermostatImplemented[3] == 1) {
         Thermostat = Langevin;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
-        if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
+        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
           cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
         } else {
@@ -324 +324 @@
       if (ThermostatImplemented[4] == 1) {
         Thermostat = Berendsen;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
       } else {
         cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
@@ -332 +332 @@
       if (ThermostatImplemented[5] == 1) {
         Thermostat = NoseHoover;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
         alpha = 0.;
       } else {
@@ -704 +704 @@
     return;
   }
+  file->close();
+  delete(file);
   RetrieveConfigPathAndName(filename);
 
@@ -722 +724 @@
   double value[3];
   
-  InitThermostats(file);
+  InitThermostats(FileBuffer);
   
   /* Namen einlesen */

src/ellipsoid.cpp

@@ -6 +6 @@
  */
 
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_vector.h>
+
+#include "boundary.hpp"
 #include "ellipsoid.hpp"
 
@@ -221 +225 @@
   set<int>::iterator current;
   int index;
-  atom *Candidate = NULL;
-  LinkedAtoms *List = NULL;
+  TesselPoint *Candidate = NULL;
+  LinkedNodes *List = NULL;
   *out << Verbose(2) << "Begin of PickRandomPointSet" << endl;
 
@@ -288 +292 @@
 //            else
 //              *out << Verbose(2) << "Cell is empty ... " << endl;
-            for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
               if ((current != PickedAtomNrs.end()) && (*current == index)) {
                 Candidate = (*Runner);
                 *out << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl;
-                x[PointsPicked++].CopyVector(&(Candidate->x));    // we have one more atom picked
+                x[PointsPicked++].CopyVector(Candidate->node);    // we have one more atom picked
                 current++;    // next pre-picked atom
               }
@@ -338 +342 @@
       //*out << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
       if (value > threshold) {
-        x[PointsPicked].CopyVector(&(Runner->second->node->x));
+        x[PointsPicked].CopyVector(Runner->second->node->node);
         PointsPicked++;
         //*out << "IN." << endl;
@@ -375 +379 @@
   Center.Zero();
   for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
-    Center.AddVector(&Runner->second->node->x);
+    Center.AddVector(Runner->second->node->node);
   Center.Scale(1./T->PointsOnBoundaryCount);
   *out << Verbose(1) << "Center is at " << Center << "." << endl;

src/ellipsoid.hpp

@@ -14 +14 @@
 #endif
 
-#include <cstdlib>
-
-#include "boundary.hpp"
+#include "linkedcell.hpp"
 #include "vector.hpp"
 

src/linkedcell.cpp

@@ -1 +1 @@
 #include "linkedcell.hpp"
 #include "molecules.hpp"
+#include "tesselation.hpp"
+
+// ========================================================= class LinkedCell ===========================================
+
 
 /** Constructor for class LinkedCell.
@@ -16 +20 @@
 
 /** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
- * \param *mol molecule structure with all Atom's
+ * \param *set LCNodeSet class with all LCNode's
  * \param RADIUS edge length of cells
  */
-LinkedCell::LinkedCell(molecule *mol, double radius)
+LinkedCell::LinkedCell(PointCloud *set, double radius)
 {
-  atom *Walker = NULL;
+  TesselPoint *Walker = NULL;
 
   RADIUS = radius;
@@ -31 +35 @@
   min.Zero();
   cout << Verbose(1) << "Begin of LinkedCell" << endl;
-  if (mol->start->next == mol->end) {
-    cerr << "ERROR: molecule contains no atoms!" << endl;
+  if (set->IsEmpty()) {
+    cerr << "ERROR: set contains no linked cell nodes!" << endl;
     return;
   }
   // 1. find max and min per axis of atoms
-  Walker = mol->start->next;
+  set->GoToFirst();
+  Walker = set->GetPoint();
   for (int i=0;i<NDIM;i++) {
-    max.x[i] = Walker->x.x[i];
-    min.x[i] = Walker->x.x[i];
+    max.x[i] = Walker->node->x[i];
+    min.x[i] = Walker->node->x[i];
   }
-  while (Walker != mol->end) {
+  set->GoToFirst();
+  while (!set->IsLast()) {
+    Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
-      if (max.x[i] < Walker->x.x[i])
-        max.x[i] = Walker->x.x[i];
-      if (min.x[i] > Walker->x.x[i])
-        min.x[i] = Walker->x.x[i];
+      if (max.x[i] < Walker->node->x[i])
+        max.x[i] = Walker->node->x[i];
+      if (min.x[i] > Walker->node->x[i])
+        min.x[i] = Walker->node->x[i];
     }
-    Walker = Walker->next;
+    set->GoToNext();
   }
   cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
 
-  // 2. find then umber of cells per axis
+  // 2. find then number of cells per axis
   for (int i=0;i<NDIM;i++) {
     N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
@@ -64 +71 @@
     return;
   }
-  LC = new LinkedAtoms[N[0]*N[1]*N[2]];
+  LC = new LinkedNodes[N[0]*N[1]*N[2]];
   for (index=0;index<N[0]*N[1]*N[2];index++) {
     LC [index].clear();
@@ -72 +79 @@
   // 4. put each atom into its respective cell
   cout << Verbose(2) << "Filling cells ... ";
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
+  set->GoToFirst();
+  while (!set->IsLast()) {
+    Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
-      n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS);
+      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
     }
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
     LC[index].push_back(Walker);
     //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
+    set->GoToNext();
   }
   cout << "done."  << endl;
@@ -118 +126 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
  */
-LinkedAtoms* LinkedCell::GetCurrentCell()
+LinkedNodes* LinkedCell::GetCurrentCell()
 {
   if (CheckBounds()) {
@@ -128 +136 @@
 };
 
-/** Calculates the index for a given atom *Walker.
- * \param *Walker atom to set index to
+/** Calculates the index for a given LCNode *Walker.
+ * \param *Walker LCNode to set index tos
  * \return if the atom is also found in this cell - true, else - false
  */
-bool LinkedCell::SetIndexToAtom(atom *Walker)
+bool LinkedCell::SetIndexToNode(TesselPoint *Walker)
 {
   bool status = false;
   for (int i=0;i<NDIM;i++) {
-    n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS);
+    n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
   }
   index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
   if (CheckBounds()) {
-    for (LinkedAtoms::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
+    for (LinkedNodes::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
       status = status || ((*Runner) == Walker);
     return status;
   } else {
-    cerr << Verbose(1) << "ERROR: Atom "<< *Walker << " at " << Walker->x << " is out of bounds." << endl;
+    cerr << Verbose(1) << "ERROR: Node at " << *Walker << " is out of bounds." << endl;
     return false;
   }

src/linkedcell.hpp

@@ +1 @@
+/*
+ * linkedcell.hpp
+ *
+ *  If the linked cell should be usable, the class has to inherit LCNodeSet and the nodes (containing the Vectors) have to inherit LCNode. This works well
+ *  for molecule and atom classes.
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
 #ifndef LINKEDCELL_HPP_
 #define LINKEDCELL_HPP_
+
+using namespace std;
 
 // include config.h
@@ -7 +19 @@
 #endif
 
-#include "molecules.hpp"
+#include <list>
 
-#define LinkedAtoms list <atom *>
+#include "defs.hpp"
+#include "helpers.hpp"
+#include "vector.hpp"
 
-class LinkedCell{
+class TesselPoint;
+class PointCloud;
+
+#define LinkedNodes list<TesselPoint *>
+
+/** Linked Cell class for containing Vectors in real space efficiently.
+ */
+class LinkedCell {
   public:
     Vector max;       // upper boundary
     Vector min;       // lower boundary
-    LinkedAtoms *LC;  // linked cell list
+    LinkedNodes *LC;  // linked cell list
     double RADIUS;    // cell edge length
     int N[NDIM];      // number of cells per axis
@@ -22 +43 @@
 
     LinkedCell();
-    LinkedCell(molecule *mol, double RADIUS);
+    LinkedCell(PointCloud *set, double RADIUS);
     ~LinkedCell();
-    LinkedAtoms* GetCurrentCell();
-    bool SetIndexToAtom(atom *Walker);
+    LinkedNodes* GetCurrentCell();
+    bool SetIndexToNode(TesselPoint *Walker);
     bool SetIndexToVector(Vector *x);
     bool CheckBounds();
 
     // not implemented yet
-    bool AddAtom(atom *Walker);
-    bool DeleteAtom(atom *Walker);
-    bool MoveAtom(atom *Walker);
+    bool AddNode(Vector *Walker);
+    bool DeleteNode(Vector *Walker);
+    bool MoveNode(Vector *Walker);
 };
 

src/molecules.cpp

@@ -42 +42 @@
   end->father = NULL;
   link(start,end);
+  InternalPointer = start;
   // init bond chain list
   first = new bond(start, end, 1, -1);
@@ -82 +83 @@
   delete(end);
   delete(start);
+};
+
+
+/** Determine center of all atoms.
+ * \param *out output stream for debugging
+ * \return pointer to allocated with central coordinates
+ */
+Vector *molecule::GetCenter(ofstream *out)
+{
+  Vector *center = DetermineCenterOfAll(out);
+  return center;
+};
+
+/** Return current atom in the list.
+ * \return pointer to atom or NULL if none present
+ */
+TesselPoint *molecule::GetPoint()
+{
+  if ((InternalPointer != start) && (InternalPointer != end))
+    return InternalPointer;
+  else
+    return NULL;
+};
+
+/** Return pointer to one after last atom in the list.
+ * \return pointer to end marker
+ */
+TesselPoint *molecule::GetTerminalPoint()
+{
+  return end;
+};
+
+/** Go to next atom.
+ * Stops at last one.
+ */
+void molecule::GoToNext()
+{
+  if (InternalPointer->next != end)
+    InternalPointer = InternalPointer->next;
+};
+
+/** Go to previous atom.
+ * Stops at first one.
+ */
+void molecule::GoToPrevious()
+{
+  if (InternalPointer->previous != start)
+    InternalPointer = InternalPointer->previous;
+};
+
+/** Goes to first atom.
+ */
+void molecule::GoToFirst()
+{
+  InternalPointer = start->next;
+};
+
+/** Goes to last atom.
+ */
+void molecule::GoToLast()
+{
+  InternalPointer = end->previous;
+};
+
+/** Checks whether we have any atoms in molecule.
+ * \return true - no atoms, false - not empty
+ */
+bool molecule::IsEmpty()
+{
+  return (start->next == end);
+};
+
+/** Checks whether we are at the last atom
+ * \return true - current atom is last one, false - is not last one
+ */
+bool molecule::IsLast()
+{
+  return (InternalPointer->next == end);
 };
 

src/molecules.hpp

@@ -25 +25 @@
 #include <vector>
 
+#include "atom.hpp"
+#include "bond.hpp"
 #include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "parser.hpp"
 #include "periodentafel.hpp"
 #include "stackclass.hpp"
+#include "tesselation.hpp"
 #include "vector.hpp"
 
-class atom;
-class bond;
 class config;
 class molecule;
@@ -38 +40 @@
 class MoleculeListClass;
 class Verbose;
-class Tesselation;
 
 /******************************** Some definitions for easier reading **********************************/
@@ -58 +59 @@
 #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 CandidateList list <class CandidateForTesselation *>
-
-#define LineMap multimap < int, class BoundaryLineSet * >
-#define LinePair pair < int, class BoundaryLineSet * >
-#define LineTestPair pair < LineMap::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> >
-
 #define MoleculeList list <molecule *>
 #define MoleculeListTest pair <MoleculeList::iterator, bool>
@@ -125 +110 @@
 };
 
-/** Single atom.
- * Class incoporates position, type
- */
-class atom {
-  public:
-    Vector x;       //!< coordinate array of atom, giving position within cell
-    Vector v;       //!< velocity array of atom
-    element *type;  //!< pointing to element
-    atom *previous; //!< previous atom in molecule list
-    atom *next;     //!< next atom in molecule list
-    atom *father;   //!< In many-body bond order fragmentations points to originating atom
-    atom *Ancestor; //!< "Father" in Depth-First-Search
-    char *Name;      //!< unique name used during many-body bond-order fragmentation
-    int FixedIon;   //!< config variable that states whether forces act on the ion or not
-    int *sort;      //!< sort criteria
-    int nr;         //!< continuous, unique number
-    int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
-    int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
-    int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
-    bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
-    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
-    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
-    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
-
-  atom();
-  atom(class atom *pointer);
-  ~atom();
-
-  bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
-  bool OutputXYZLine(ofstream *out) const;
-  atom *GetTrueFather();
-  bool Compare(atom &ptr);
-
-  private:
-};
-
-ostream & operator << (ostream &ost, const atom &a);
-
-/** Bonds between atoms.
- * Class incorporates bonds between atoms in a molecule,
- * used to derive tge fragments in many-body bond order
- * calculations.
- */
-class bond {
-  public:
-    atom *leftatom;    //!< first bond partner
-    atom *rightatom;  //!< second bond partner
-    bond *previous; //!< previous atom in molecule list
-    bond *next;     //!< next atom in molecule list
-    int HydrogenBond;  //!< Number of hydrogen atoms in the bond
-    int BondDegree;    //!< single, double, triple, ... bond
-    int nr;           //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
-    bool Cyclic;      //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
-    enum EdgeType Type;//!< whether this is a tree or back edge
-
-  atom * GetOtherAtom(atom *Atom) const;
-  bond * GetFirstBond();
-  bond * GetLastBond();
-
-  bool MarkUsed(enum Shading color);
-  enum Shading IsUsed();
-  void ResetUsed();
-  bool Contains(const atom *ptr);
-  bool Contains(const int nr);
-
-  bond();
-  bond(atom *left, atom *right);
-  bond(atom *left, atom *right, int degree);
-  bond(atom *left, atom *right, int degree, int number);
-  ~bond();
-
-  private:
-    enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
-};
-
-
-ostream & operator << (ostream &ost, const bond &b);
 
 #define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
@@ -210 +118 @@
  * Class incorporates number of types
  */
-class molecule {
+class molecule : public PointCloud {
   public:
     double cell_size[6];//!< cell size
@@ -234 +142 @@
     char name[MAXSTRINGSIZE];         //!< arbitrary name
     int IndexNr;        //!< index of molecule in a MoleculeListClass
+    class Tesselation *TesselStruct;
 
   molecule(periodentafel *teil);
   ~molecule();
+
+  // re-definition of virtual functions from PointCloud
+  Vector *GetCenter(ofstream *out);
+  TesselPoint *GetPoint();
+  TesselPoint *GetTerminalPoint();
+  void GoToNext();
+  void GoToPrevious();
+  void GoToFirst();
+  void GoToLast();
+  bool IsEmpty();
+  bool IsLast();
 
   /// remove atoms from molecule.
@@ -350 +270 @@
   private:
   int last_atom;      //!< number given to last atom
+  atom *InternalPointer;  //!< internal pointer for PointCloud
 };
 
@@ -524 +445 @@
   char *GetDefaultPath() const;
   void SetDefaultPath(const char *path);
-  void InitThermostats(ifstream *source);
+  void InitThermostats(class ConfigFileBuffer *fb);
 };
 

src/vector.hpp

@@ -1 +1 @@
 #ifndef VECTOR_HPP_
 #define VECTOR_HPP_
+
+using namespace std;
+
+#include "helpers.hpp"
 
 class Vector;