Index: src/#border.cpp#
===================================================================
--- src/#border.cpp#	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
+++ src/#border.cpp#	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -0,0 +1,302 @@
+#include "molecules.hpp"
+#include "boundary.hpp"
+
+
+void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol)
+{
+  /* d2 ist der Normalenvektor auf dem Dreieck,
+   * d1 ist der Vektor, der normal auf der Kante und d2 steht.
+   */
+  Vector dif_a;
+  Vector dif_b;
+  Vector Chord;
+  Vector AngleCheck;
+  atom *Walker;
+
+  dif_a.CopyVector(&a.x);
+  dif_a.SubtractVector(&Candidate.x);
+  dif_b.CopyVector(&b.x);
+  dif_b.SubtractVector(&Candidate.x);
+  Chord.CopyVector(&a.x);
+  Chord.SubtractVector(&b.x);
+  AngleCheck.CopyVector(&dif_a);
+  AngleCheck.ProjectOntoPlane(&Chord);
+
+  //Storage eintrag fuer aktuelles Atom
+  if (Chord.Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  {
+
+    if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
+      {
+	Storage[0]=(double)Candidate.nr;
+	Storage[1]=1;
+	Storage[2]=AngleCheck.Angle(d2);
+      }
+    else
+      {
+	if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
+	    (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
+	  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+	  {
+	    Storage[0]=(double)Candidate.nr;
+	    Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+	    Storage[2]=AngleCheck.Angle(d2);
+	  }
+      }
+  }
+
+  if (n<5)
+    {
+      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
+	{
+	  while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
+	    {
+	      Walker = Walker->next;
+	    }
+
+	  Find_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
+	}
+    }
+};
+
+
+void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
+{
+  Vector CenterOfLine = Line.endpoints[0]->node->x;
+  Vector direction1;
+  Vector direction2;
+  Vector helper;
+  atom* Walker;
+
+  double Storage[3];
+  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
+  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
+
+  
+  helper.CopyVector(&(Line.endpoints[0]->node->x));
+  for (int i =0; i<3; i++)
+    {
+      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr!=Line.endpoints[1]->node->nr)
+	{
+	  helper.SubtractVector(&T.endpoints[i]->node->x);
+	  break;
+	}
+    }
+
+
+  direction1.CopyVector(&Line.endpoints[0]->node->x);
+  direction1.SubtractVector(&Line.endpoints[1]->node->x);
+  direction1.VectorProduct(T.NormalVector);
+
+  if (direction1.ScalarProduct(&helper)>0)
+    {
+      direction1.Scale(-1);
+    }
+
+  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &direction1, T.NormalVector, Storage, RADIUS, *mol);
+  
+  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
+  // Next Triangle is Line, atom with number in Storage[0]
+
+  Walker= mol->start;
+  while (Walker->nr != (int)Storage[0])
+    {
+      Walker = Walker->next;
+    }
+
+  AddPoint(Walker);
+
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BLS[2] = new class BoundaryLineSet(Line);
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
+    {
+      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
+	{
+	  LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+	  LinesOnBoundaryCount++;
+	}
+    }
+  BTS->GetNormalVector(*BTS->NormalVector);
+
+  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
+      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
+    {
+      BTS->NormalVector->Scale(-1);
+    }
+  
+};
+
+
+void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
+{
+  int i;
+  Vector *AngleCheck;
+  atom* Walker;
+
+  AngleCheck->CopyVector(&Candidate.x);
+  AngleCheck->SubtractVector(&a.x);
+  if (AngleCheck->ScalarProduct(&Oben) < Storage[1])
+    {
+      Storage[0]=(double)(Candidate.nr);
+      Storage[1]=AngleCheck->ScalarProduct(&Oben);
+    };
+
+  if (n<5)
+    {
+      for (i = 0; i< mol.NumberOfBondsPerAtom[Candidate.nr]; i++)
+	{
+	  Walker = mol.start;
+	  while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
+	    {
+	      Walker = Walker->next;
+	    };
+	  
+	  Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
+	    };
+    };
+  
+
+};
+
+
+void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
+{
+  int i=0;
+  atom Walker;
+  atom Walker2;
+  atom Walker3;
+  int max_index[3];
+  double max_coordinate[3];
+  Vector Oben;
+  Vector helper;
+
+  Oben.Zero();
+
+
+  for(i =0; i<3; i++)
+    {
+      max_index[i] =-1;
+      max_coordinate[i] =-1;
+    }
+
+  Walker = *mol.start;
+  while (Walker.next != NULL)
+    {
+      for (i=0; i<3; i++)
+	{
+	  if (Walker.x.x[i]>max_coordinate[i])
+	    {
+	      max_coordinate[i]=Walker.x.x[i];
+	      max_index[i]=Walker.nr;
+	    }
+	}
+    }
+
+  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
+  const int k=0;
+
+  Oben.x[k]=1;
+  Walker = *mol.start;
+  while (Walker.nr != max_index[k])
+    {
+      Walker = *Walker.next;
+    }
+
+  double Storage[3];
+  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
+  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[2]=-10.;  // This will be an angle looking for the third point.
+
+
+  for (i=0; i< mol.NumberOfBondsPerAtom[Walker.nr]; i++)
+    {
+      Walker2 = *mol.start;
+      while (Walker2.nr != (mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom->nr : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr) ) 
+	{
+	  Walker2 = *Walker2.next;
+	}
+
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
+    }
+
+  Walker2 = *mol.start;
+
+  while (Walker2.nr != int(Storage[0]))
+    {
+      Walker = *Walker.next;
+    }
+  
+  helper.CopyVector(&Walker.x);
+  helper.SubtractVector(&Walker2.x);
+  Oben.ProjectOntoPlane(&helper);
+  helper.VectorProduct(&Oben);
+
+       Find_next_suitable_point(Walker, Walker2, *(mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom), 0, &helper, &Oben, Storage, RADIUS, mol);
+  Walker3 = *mol.start;
+  while (Walker3.nr != int(Storage[0]))
+    {
+      Walker3 = *Walker3.next;
+    }
+
+  //Starting Triangle is Walker, Walker2, index Storage[0]
+
+  AddPoint(&Walker);
+  AddPoint(&Walker2);
+  AddPoint(&Walker3);
+
+  BPS[0] = new class BoundaryPointSet(&Walker);
+  BPS[1] = new class BoundaryPointSet(&Walker2);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(&Walker);
+  BPS[1] = new class BoundaryPointSet(&Walker3);
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(&Walker);
+  BPS[1] = new class BoundaryPointSet(&Walker2);
+  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
+
+  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++;
+    };
+
+       BTS->GetNormalVector(*BTS->NormalVector);
+
+       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
+	 {
+	   BTS->NormalVector->Scale(-1);
+	 }
+};
+
+
+void Find_non_convex_border(Tesselation* Tess, molecule mol) 
+{
+  const double RADIUS =6;
+  Tess->Find_starting_triangle(mol, RADIUS);
+
+  for (LineMap::iterator baseline = Tess->LinesOnBoundary.begin(); baseline != Tess->LinesOnBoundary.end(); baseline++) 
+    if (baseline->second->TrianglesCount == 1) 
+      {
+	Tess->Find_next_suitable_triangle(&mol, *(baseline->second), baseline->second->triangles.begin()->second, RADIUS); //the line is there, so there is a triangle, but only one.
+
+      }
+    else
+      {
+	printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+      }
+};
Index: src/#makefile#
===================================================================
--- src/#makefile#	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
+++ src/#makefile#	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -0,0 +1,10 @@
+CC = gcc
+CFLAGS = -g -Wall
+CompileFlags = -lgsl -lgslcblas -lm
+
+
+.c.o:
+	$(CC) $(CFLAGS) -c $<
+
+clean:
+	rm -f *.o
Index: src/#molecules.hpp#
===================================================================
--- src/#molecules.hpp#	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
+++ src/#molecules.hpp#	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -0,0 +1,463 @@
+/** \file molecules.hpp
+ *
+ * Class definitions of atom and molecule, element and periodentafel 
+ */
+
+#ifndef MOLECULES_HPP_
+#define MOLECULES_HPP_
+
+using namespace std;
+
+// GSL headers
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_heapsort.h>
+
+// STL headers
+#include <map>
+#include <set>
+#include <deque>
+#include <list>
+#include <vector>
+
+#include "helpers.hpp"
+#include "parser.hpp"
+#include "periodentafel.hpp"
+#include "stackclass.hpp"
+#include "vector.hpp"
+
+class atom;
+class bond;
+class config;
+class molecule;
+class MoleculeListClass;
+class Verbose;
+
+/******************************** Some definitions for easier reading **********************************/
+
+#define KeyStack deque<int>
+#define KeySet set<int>
+#define NumberValuePair pair<int, double>
+#define Graph map <KeySet, NumberValuePair, KeyCompare >
+#define GraphPair pair <KeySet, NumberValuePair >
+#define KeySetTestPair pair<KeySet::iterator, bool>
+#define GraphTestPair pair<Graph::iterator, bool>
+
+#define DistancePair pair < double, atom* >
+#define DistanceMap multimap < double, atom* >
+#define DistanceTestPair pair < DistanceMap::iterator, bool>
+
+#define Boundaries map <double, DistancePair >
+#define BoundariesPair pair<double, DistancePair >
+#define BoundariesTestPair pair< Boundaries::iterator, bool>
+
+#define PointMap map < int, class BoundaryPointSet * > 
+#define PointPair pair < int, class BoundaryPointSet * > 
+#define PointTestPair pair < PointMap::iterator, bool > 
+
+#define LineMap map < int, class BoundaryLineSet * > 
+#define LinePair pair < int, class BoundaryLineSet * > 
+#define LineTestPair pair < LinePair::iterator, bool > 
+
+#define TriangleMap map < int, class BoundaryTriangleSet * > 
+#define TrianglePair pair < int, class BoundaryTriangleSet * > 
+#define TriangleTestPair pair < TrianglePair::iterator, bool > 
+
+#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
+#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
+
+/******************************** Some small functions and/or structures **********************************/
+
+struct KeyCompare
+{
+  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
+};
+
+struct Trajectory
+{
+  vector<Vector> R;  //!< position vector
+  vector<Vector> U;  //!< velocity vector
+  vector<Vector> F;  //!< last force vector
+  atom *ptr;         //!< pointer to atom whose trajectory we contain
+};
+
+//bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)
+inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
+inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph 
+int CompareDoubles (const void * a, const void * b);
+
+
+/************************************* Class definitions ****************************************/
+
+
+// some algebraic matrix stuff
+#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])  //!< hard-coded determinant of a 3x3 matrix
+#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix
+
+
+/** Parameter structure for least square minimsation.
+ */
+struct LSQ_params {
+  Vector **vectors;
+  int num;
+};
+
+double LSQ(const gsl_vector * x, void * params);
+
+/** Parameter structure for least square minimsation.
+ */
+struct lsq_params {
+  gsl_vector *x;
+  const molecule *mol;
+  element *type;
+};
+
+
+
+/** 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();
+  
+  bool Output(int ElementNo, int AtomNo, ofstream *out) const;
+  bool OutputXYZLine(ofstream *out) const;
+  atom *GetTrueFather();
+  bool Compare(atom &ptr);
+  
+  private:
+};
+
+ostream & operator << (ostream &ost, 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, bond &b);
+
+class MoleculeLeafClass;
+
+/** The complete molecule.
+ * Class incorporates number of types
+ */
+class molecule {
+  public:
+    double cell_size[6];//!< cell size
+    periodentafel *elemente; //!< periodic table with each element
+    atom *start;        //!< start of atom list
+    atom *end;          //!< end of atom list
+    bond *first;        //!< start of bond list
+    bond *last;         //!< end of bond list
+    bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
+    map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
+    int MDSteps;        //!< The number of MD steps in Trajectories
+    int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
+    int AtomCount;					//!< number of atoms, brought up-to-date by CountAtoms()
+    int BondCount;					//!< number of atoms, brought up-to-date by CountBonds()
+    int ElementCount;       //!< how many unique elements are therein
+    int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
+    int NoNonHydrogen;	//!< number of non-hydrogen atoms in molecule
+    int NoNonBonds;     //!< number of non-hydrogen bonds in molecule
+    int NoCyclicBonds;  //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
+    double BondDistance;  //!< typical bond distance used in CreateAdjacencyList() and furtheron
+  
+  molecule(periodentafel *teil);
+  ~molecule();
+  
+  /// remove atoms from molecule.
+  bool AddAtom(atom *pointer);
+  bool RemoveAtom(atom *pointer);
+  bool CleanupMolecule();
+
+  /// Add/remove atoms to/from molecule.
+  atom * AddCopyAtom(atom *pointer);
+  bool AddXYZFile(string filename);
+  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);  
+  bond * AddBond(atom *first, atom *second, int degree);
+  bool RemoveBond(bond *pointer);
+  bool RemoveBonds(atom *BondPartner);
+    
+  /// Find atoms.
+  atom * FindAtom(int Nr) const; 
+  atom * AskAtom(string text);
+
+  /// Count and change present atoms' coordination.
+  void CountAtoms(ofstream *out);
+  void CountElements();
+  void CalculateOrbitals(class config &configuration);
+  bool CenterInBox(ofstream *out, Vector *BoxLengths);
+  void CenterEdge(ofstream *out, Vector *max); 
+  void CenterOrigin(ofstream *out, Vector *max); 
+  void CenterGravity(ofstream *out, Vector *max);
+  void Translate(const Vector *x);
+  void Mirror(const Vector *x);
+  void Align(Vector *n);
+  void Scale(double **factor);
+  void DetermineCenter(Vector &center);
+  Vector * DetermineCenterOfGravity(ofstream *out);
+  Vector * DetermineCenterOfAll(ofstream *out);
+  void SetBoxDimension(Vector *dim);
+  double * ReturnFullMatrixforSymmetric(double *cell_size);
+  void ScanForPeriodicCorrection(ofstream *out);
+	void PrincipalAxisSystem(ofstream *out, bool DoRotate);
+	double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
+	bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
+	
+  bool CheckBounds(const Vector *x) const;
+  void GetAlignvector(struct lsq_params * par) const;
+
+  /// Initialising routines in fragmentation  
+  void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
+  void CreateListOfBondsPerAtom(ofstream *out);
+  
+  // Graph analysis
+  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
+  void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
+  bool PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack);
+  bond * FindNextUnused(atom *vertex);
+  void SetNextComponentNumber(atom *vertex, int nr);
+  void InitComponentNumbers();
+  void OutputComponentNumber(ofstream *out, atom *vertex);
+  void ResetAllBondsToUnused();
+  void ResetAllAtomNumbers();
+  int CountCyclicBonds(ofstream *out);
+  bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
+  string GetColor(enum Shading color);
+
+  molecule *CopyMolecule();
+  
+  /// Fragment molecule by two different approaches:
+  int FragmentMolecule(ofstream *out, int Order, config *configuration);
+  bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
+  bool StoreAdjacencyToFile(ofstream *out, char *path);
+  bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
+  bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
+  bool StoreOrderAtSiteFile(ofstream *out, char *path);
+  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
+  bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
+  bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
+  bool CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex);
+  bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
+  void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
+  /// -# BOSSANOVA
+  void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
+	int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
+  bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
+  molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
+  void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
+  int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
+  int GuesstimateFragmentCount(ofstream *out, int order);
+  	  
+  // Recognize doubly appearing molecules in a list of them   
+  int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
+  int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
+	
+  // Output routines.
+  bool Output(ofstream *out);
+  bool OutputTrajectories(ofstream *out);
+  void OutputListOfBonds(ofstream *out) const;
+  bool OutputXYZ(ofstream *out) const;
+  bool OutputTrajectoriesXYZ(ofstream *out);
+  bool Checkout(ofstream *out) const;
+  bool OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output);
+
+  private:
+  int last_atom;      //!< number given to last atom
+};
+
+/** A list of \a molecule classes.
+ */
+class MoleculeListClass {
+  public:
+    molecule **ListOfMolecules;   //!< pointer list of fragment molecules to check for equality
+    int NumberOfMolecules;        //!< Number of entries in \a **FragmentList and of to be returned one.
+    int NumberOfTopAtoms;         //!< Number of atoms in the molecule from which all fragments originate
+    
+  MoleculeListClass();
+  MoleculeListClass(int Num, int NumAtoms);
+  ~MoleculeListClass();
+
+  /// Output configs.
+  bool AddHydrogenCorrection(ofstream *out, char *path);
+  bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
+  bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
+  void Output(ofstream *out);
+  
+  private:
+};
+
+
+/** A leaf for a tree of \a molecule class
+ * Wraps molecules in a tree structure
+ */
+class MoleculeLeafClass {
+  public:
+    molecule *Leaf;                   //!< molecule of this leaf 
+    //MoleculeLeafClass *UpLeaf;        //!< Leaf one level up
+    //MoleculeLeafClass *DownLeaf;      //!< First leaf one level down
+    MoleculeLeafClass *previous;  //!< Previous leaf on this level
+    MoleculeLeafClass *next;      //!< Next leaf on this level
+
+  //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
+  MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
+  ~MoleculeLeafClass();
+
+  bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
+  bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
+  bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
+  bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
+  bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
+  void TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
+  int Count() const;
+};
+
+/** The config file.
+ * The class contains all parameters that control a dft run also functions to load and save.
+ */
+class config {
+  public:
+    int PsiType;
+    int MaxPsiDouble;
+    int PsiMaxNoUp;
+    int PsiMaxNoDown;
+    int MaxMinStopStep;
+    int InitMaxMinStopStep;
+    int ProcPEGamma;
+    int ProcPEPsi;
+    char *configpath;
+    char *configname;
+    bool FastParsing;
+    double Deltat;
+    
+    private:
+    char *mainname;
+    char *defaultpath;
+    char *pseudopotpath;
+    
+    int DoOutVis;
+    int DoOutMes;
+    int DoOutNICS;
+    int DoOutOrbitals;
+    int DoOutCurrent;
+    int DoFullCurrent;
+    int DoPerturbation;
+    int DoWannier;
+    int CommonWannier;
+    double SawtoothStart;
+    int VectorPlane;
+    double VectorCut;
+    int UseAddGramSch;
+    int Seed;
+    
+    int MaxOuterStep;
+    int OutVisStep;
+    int OutSrcStep;
+    double TargetTemp;
+    int ScaleTempStep;
+    int MaxPsiStep;
+    double EpsWannier;
+    
+    int MaxMinStep;
+    double RelEpsTotalEnergy;
+    double RelEpsKineticEnergy;
+    int MaxMinGapStopStep;
+    int MaxInitMinStep;
+    double InitRelEpsTotalEnergy;
+    double InitRelEpsKineticEnergy;
+    int InitMaxMinGapStopStep;
+    
+    //double BoxLength[NDIM*NDIM];
+    
+    double ECut;
+    int MaxLevel;
+    int RiemannTensor;
+    int LevRFactor;
+    int RiemannLevel;
+    int Lev0Factor;
+    int RTActualUse;
+    int AddPsis;
+    
+    double RCut;
+    int StructOpt;
+    int IsAngstroem;
+    int RelativeCoord;
+    int MaxTypes;
+
+  
+  int ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical);
+  
+  public:
+  config();
+  ~config();
+
+  int TestSyntax(char *filename, periodentafel *periode, molecule *mol);
+  void Load(char *filename, periodentafel *periode, molecule *mol);
+  void LoadOld(char *filename, periodentafel *periode, molecule *mol);
+  void RetrieveConfigPathAndName(string filename);
+  bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
+  bool SaveMPQC(const char *filename, molecule *mol) const;
+  void Edit(molecule *mol);
+  bool GetIsAngstroem() const;
+  char *GetDefaultPath() const;
+  void SetDefaultPath(const char *path);
+};
+
+#endif /*MOLECULES_HPP_*/
+
Index: src/.#border.cpp
===================================================================
--- src/.#border.cpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
+++ src/.#border.cpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -0,0 +1,1 @@
+neuen@stud64dc-06.wiss-stud.31428:1225702725
Index: src/.#molecules.hpp
===================================================================
--- src/.#molecules.hpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
+++ src/.#molecules.hpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -0,0 +1,1 @@
+neuen@stud64dc-06.wiss-stud.31192:1225702725
Index: src/TAGS
===================================================================
--- src/TAGS	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
+++ src/TAGS	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -0,0 +1,376 @@
+
+
+analyzer.cpp,17
+int main(23,516
+
+
+atom.cpp,262
+atom::atom(13,258
+atom::~atom(35,653
+atom *atom::GetTrueFather(atom::GetTrueFather45,973
+bool atom::Output(61,1416
+bool atom::OutputXYZLine(78,2040
+ostream & operator << operator <<87,2262
+bool atom::Compare(97,2547
+bool operator < operator <105,2652
+
+
+bond.cpp,385
+bond::bond(14,290
+bond::bond(34,691
+bond::bond(52,1136
+bond::~bond(73,1597
+ostream & operator << operator <<84,1788
+atom * bond::GetOtherAtom(94,2199
+bond * bond::GetFirstBond(107,2583
+bond * bond::GetLastBond(116,2865
+enum Shading bond::IsUsed(124,3000
+bool bond::Contains(133,3217
+bool bond::Contains(142,3475
+bool bond::MarkUsed(150,3691
+void bond::ResetUsed(163,3973
+
+
+border.cpp,229
+void Find_next_suitable_point(5,51
+void Tesselation::Find_next_suitable_triangle(63,2078
+void Find_second_point_for_Tesselation(141,4691
+void Tesselation::Find_starting_triangle(173,5558
+void Find_non_convex_border(287,8778
+
+
+boundary.cpp,1347
+BoundaryPointSet::BoundaryPointSet(12,397
+BoundaryPointSet::BoundaryPointSet(18,469
+BoundaryPointSet::~BoundaryPointSet(25,578
+void BoundaryPointSet::AddLine(31,704
+ostream & operator << operator <<42,1045
+BoundaryLineSet::BoundaryLineSet(50,1273
+BoundaryLineSet::BoundaryLineSet(58,1396
+BoundaryLineSet::~BoundaryLineSet(72,1836
+void BoundaryLineSet::AddTriangle(87,2407
+ostream & operator << operator <<94,2656
+BoundaryTriangleSet::BoundaryTriangleSet(103,2939
+BoundaryTriangleSet::BoundaryTriangleSet(112,3075
+BoundaryTriangleSet::~BoundaryTriangleSet(144,4167
+void BoundaryTriangleSet::GetNormalVector(159,4691
+ostream & operator << operator <<169,5082
+class BoundaryPointSet * GetCommonEndpoint(182,5559
+Boundaries * GetBoundaryPoints(209,6893
+double * GetDiametersOfCluster(384,15165
+double VolumeOfConvexEnvelope(470,20313
+void PrepareClustersinWater(610,27033
+Tesselation::Tesselation(692,31276
+Tesselation::~Tesselation(702,31494
+void Tesselation::GuessStartingTriangle(714,31995
+void Tesselation::TesselateOnBoundary(849,39116
+void Tesselation::AddPoint(1003,49558
+void Find_next_suitable_point(1024,50156
+void Tesselation::Find_next_suitable_triangle(1082,52183
+void Find_second_point_for_Tesselation(1160,54796
+void Tesselation::Find_starting_triangle(1192,55663
+void Find_non_convex_border(1306,58883
+
+
+builder.cpp,378
+static void AddAtoms(62,2281
+static void CenterAtoms(259,10538
+static void AlignAtoms(314,12598
+static void MirrorAtoms(382,15117
+static void RemoveAtoms(429,16855
+static void MeasureAtoms(487,18852
+static void FragmentAtoms(598,23710
+static void testroutine(619,24543
+static void SaveConfig(692,27108
+static int ParseCommandLineOptions(771,29725
+int main(1303,56039
+
+
+config.cpp,445
+config::config(13,261
+config::~config(89,2039
+void config::Edit(102,2590
+int config::TestSyntax(367,16198
+bool config::GetIsAngstroem(389,16806
+char * config::GetDefaultPath(397,16955
+void config::SetDefaultPath(406,17101
+void config::RetrieveConfigPathAndName(414,17286
+void config::Load(444,18309
+void config::LoadOld(773,36012
+bool config::Save(977,45063
+bool config::SaveMPQC(1084,52722
+int config::ParseForParameter(1213,58282
+
+
+datacreator.cpp,954
+bool OpenOutputFile(19,468
+bool AppendOutputFile(37,1009
+bool CreateDataEnergyOrder(56,1690
+bool CreateDataDeltaEnergyOrder(89,3244
+bool CreateDataForcesOrder(126,5185
+bool CreateDataDeltaForcesOrder(158,6656
+bool CreateDataDeltaForcesOrderPerAtom(190,8202
+bool CreateDataForcesOrderPerAtom(237,10288
+bool CreateDataFragment(266,11393
+void CreateMaxFragmentOrder(294,12696
+void CreateMinFragmentOrder(311,13550
+bool CreateDataFragmentOrder(331,14537
+void CreateEnergy(358,15708
+void CreateMinimumForce(369,16220
+void CreateMeanForce(403,17598
+void CreateMaximumForce(428,18520
+void CreateSameForce(450,19299
+void CreateVectorSumForce(460,19630
+void CreatePlotHeader(480,20389
+bool CreatePlotOrder(517,21999
+void AbsEnergyPlotLine(538,23154
+void EnergyPlotLine(562,23997
+void ForceMagnitudePlotLine(586,24845
+void AbsFirstForceValuePlotLine(617,26051
+void BoxesForcePlotLine(648,27192
+void BoxesFirstForceValuePlotLine(680,28575
+
+
+element.cpp,112
+element::element(13,270
+element::~element(23,412
+bool element::Output(28,507
+bool element::Checkout(43,1115
+
+
+helpers.cpp,325
+double ask_value(15,316
+void debug_in(30,757
+void debug_in(34,918
+void * Malloc(42,1223
+void * Calloc(56,1671
+void * ReAlloc(71,2182
+void Free(88,2787
+char* MallocString(103,3192
+void bound(121,3685
+void flip(134,3996
+int pot(147,4205
+char *FixedDigitNumber(FixedDigitNumber161,4624
+bool IsValidNumber(188,5558
+
+
+joiner.cpp,17
+int main(17,429
+
+
+moleculelist.cpp,859
+MoleculeListClass::MoleculeListClass(13,290
+MoleculeListClass::MoleculeListClass(21,496
+MoleculeListClass::~MoleculeListClass(33,889
+int MolCompare(52,1584
+void MoleculeListClass::Output(127,3994
+bool MoleculeListClass::AddHydrogenCorrection(142,4664
+bool MoleculeListClass::StoreForcesFile(329,12009
+bool MoleculeListClass::OutputConfigForListOfFragments(382,14293
+MoleculeLeafClass::MoleculeLeafClass(492,18874
+MoleculeLeafClass::~MoleculeLeafClass(512,19382
+bool MoleculeLeafClass::AddLeaf(550,20693
+bool MoleculeLeafClass::FillBondStructureFromReference(564,21475
+bool MoleculeLeafClass::FillRootStackForSubgraphs(622,24472
+bool MoleculeLeafClass::FillListOfLocalAtoms(662,26126
+bool MoleculeLeafClass::AssignKeySetsToFragment(696,27829
+void MoleculeLeafClass::TranslateIndicesToGlobalIDs(762,31354
+int MoleculeLeafClass::Count(786,32584
+
+
+molecules.cpp,4521
+double LSQ 16,351
+molecule::molecule(36,961
+molecule::~molecule(69,1815
+bool molecule::AddAtom(88,2550
+atom * molecule::AddCopyAtom(116,3584
+bool molecule::AddHydrogenReplacementAtom(174,7697
+bool molecule::AddXYZFile(438,21482
+molecule *molecule::CopyMolecule(molecule::CopyMolecule499,23453
+bond * molecule::AddBond(563,25434
+bool molecule::RemoveBond(582,26231
+bool molecule::RemoveBonds(594,26664
+void molecule::SetBoxDimension(603,26931
+bool molecule::CenterInBox(617,27297
+void molecule::CenterEdge(666,28615
+void molecule::CenterOrigin(702,29700
+Vector * molecule::DetermineCenterOfAll(725,30355
+Vector * molecule::DetermineCenterOfGravity(753,31098
+void molecule::CenterGravity(782,31922
+void molecule::Scale(796,32208
+void molecule::Translate(811,32523
+void molecule::Mirror(826,32863
+void molecule::DetermineCenter(841,33204
+void molecule::PrincipalAxisSystem(909,35771
+bool molecule::VerletForceIntegration(1020,39893
+void molecule::Align(1131,44421
+bool molecule::RemoveAtom(1193,46676
+bool molecule::CleanupMolecule(1208,47333
+atom * molecule::FindAtom(1217,47566
+atom * molecule::AskAtom(1231,47958
+bool molecule::CheckBounds(1250,48510
+double LeastSquareDistance 1267,49011
+void molecule::GetAlignvector(1300,50133
+bool molecule::Output(1374,52004
+bool molecule::OutputTrajectories(1409,53012
+void molecule::OutputListOfBonds(1457,55095
+bool molecule::Checkout(1480,55797
+bool molecule::OutputTrajectoriesXYZ(1488,56002
+bool molecule::OutputXYZ(1517,56942
+void molecule::CountAtoms(1544,57672
+void molecule::CountElements(1579,58991
+int molecule::CountCyclicBonds(1601,59545
+string molecule::GetColor(1630,60399
+void molecule::CalculateOrbitals(1655,60877
+void molecule::CreateAdjacencyList(1696,63214
+MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(1884,73778
+void molecule::CyclicStructureAnalysis(2078,84349
+void molecule::SetNextComponentNumber(2284,95178
+void molecule::OutputComponentNumber(2305,95927
+void molecule::InitComponentNumbers(2313,96194
+bond * molecule::FindNextUnused(2330,96839
+void molecule::ResetAllBondsToUnused(2341,97177
+void molecule::ResetAllAtomNumbers(2352,97396
+void OutputAlreadyVisited(2365,97689
+int molecule::GuesstimateFragmentCount(2382,98216
+bool molecule::ScanBufferIntoKeySet(2403,99047
+bool molecule::ParseKeySetFile(2435,100271
+bool molecule::StoreKeySetFile(2515,103341
+bool molecule::StoreAdjacencyToFile(2571,105077
+bool molecule::CheckAdjacencyFileAgainstMolecule(2606,106282
+bool molecule::CheckOrderAtSite(2684,109955
+bool molecule::CreateMappingLabelsToConfigSequence(2847,118616
+int molecule::FragmentMolecule(2892,120867
+bool molecule::PickLocalBackEdges(3084,130811
+bool molecule::StoreOrderAtSiteFile(3121,132491
+bool molecule::ParseOrderAtSiteFromFile(3152,133743
+void molecule::CreateListOfBondsPerAtom(3209,136178
+void molecule::BreadthFirstSearchAdd(3280,139568
+bool molecule::BuildInducedSubgraph(3404,147541
+int molecule::LookForRemovalCandidate(3459,149968
+molecule * molecule::StoreFragmentFromKeySet(3487,151103
+struct UniqueFragments 3782,166770
+  config *configuration;UniqueFragments::configuration3783,166795
+  atom *Root;UniqueFragments::Root3784,166820
+  Graph *Leaflet;UniqueFragments::Leaflet3785,166834
+  KeySet *FragmentSet;UniqueFragments::FragmentSet3786,166852
+  int ANOVAOrder;UniqueFragments::ANOVAOrder3787,166875
+  int FragmentCounter;UniqueFragments::FragmentCounter3788,166893
+  int CurrentIndex;UniqueFragments::CurrentIndex3789,166916
+  double TEFactor;UniqueFragments::TEFactor3790,166936
+  int *ShortestPathList;UniqueFragments::ShortestPathList3791,166955
+  bool **UsedList;UniqueFragments::UsedList3792,166980
+  bond **BondsPerSPList;UniqueFragments::BondsPerSPList3793,166999
+  int *BondsPerSPCount;UniqueFragments::BondsPerSPCount3794,167024
+void molecule::SPFragmentGenerator(3811,168011
+bool molecule::CheckForConnectedSubgraph(3945,176059
+int molecule::PowerSetGenerator(4004,178457
+void molecule::ScanForPeriodicCorrection(4156,186870
+double * molecule::ReturnFullMatrixforSymmetric(4242,190359
+bool KeyCompare::operator() KeyCompare::operator()4257,190766
+inline void InsertFragmentIntoGraph(4294,192008
+inline void InsertGraphIntoGraph(4323,193480
+void molecule::FragmentBOSSANOVA(4357,195580
+inline int CompareDoubles 4547,207047
+int * molecule::IsEqualToWithinThreshold(4563,207650
+int * molecule::GetFatherSonAtomicMap(4688,213580
+bool molecule::OutputTemperatureFromTrajectories(4734,215516
+
+
+parser.cpp,1173
+bool FilePresent(24,588
+bool TestParams(43,1017
+MatrixContainer::MatrixContainer(56,1344
+MatrixContainer::~MatrixContainer(69,1863
+bool MatrixContainer::ParseMatrix(109,3580
+bool MatrixContainer::ParseFragmentMatrix(210,7695
+bool MatrixContainer::AllocateMatrix(260,9629
+bool MatrixContainer::ResetMatrix(284,10737
+double MatrixContainer::FindMaxValue(296,11067
+double MatrixContainer::FindMinValue(312,11530
+bool MatrixContainer::SetLastMatrix(330,12035
+bool MatrixContainer::SetLastMatrix(343,12489
+bool MatrixContainer::SumSubManyBodyTerms(358,13130
+bool MatrixContainer::WriteTotalFragments(412,16895
+bool MatrixContainer::WriteLastMatrix(445,17902
+bool EnergyMatrix::ParseIndices(473,18785
+bool EnergyMatrix::SumSubEnergy(494,19731
+bool EnergyMatrix::ParseFragmentMatrix(518,20903
+bool ForceMatrix::ParseIndices(550,22485
+bool ForceMatrix::SumSubForces(597,24393
+bool ForceMatrix::ParseFragmentMatrix(628,25589
+KeySetsContainer::KeySetsContainer(679,27595
+KeySetsContainer::~KeySetsContainer(690,27804
+bool KeySetsContainer::ParseKeySets(707,28606
+bool KeySetsContainer::ParseManyBodyTerms(751,30268
+bool KeySetsContainer::Contains(814,32310
+
+
+periodentafel.cpp,516
+periodentafel::periodentafel(16,396
+periodentafel::~periodentafel(29,703
+bool periodentafel::AddElement(40,947
+bool periodentafel::RemoveElement(52,1308
+bool periodentafel::CleanupPeriodtable(60,1520
+element * periodentafel::FindElement(70,1786
+element * periodentafel::FindElement(94,2687
+element * periodentafel::AskElement(107,3041
+bool periodentafel::Output(123,3370
+bool periodentafel::Checkout(141,3816
+bool periodentafel::LoadPeriodentafel(166,4629
+bool periodentafel::StorePeriodentafel(305,9306
+
+
+vector.cpp,1457
+Vector::Vector(14,267
+Vector::Vector(18,356
+Vector::~Vector(22,484
+double Vector::Distance(28,639
+double Vector::PeriodicDistance(41,1061
+void Vector::KeepPeriodic(78,2451
+double Vector::ScalarProduct(109,3471
+void Vector::VectorProduct(124,3904
+void Vector::ProjectOntoPlane(139,4272
+double Vector::Projection(151,4570
+double Vector::Norm(159,4723
+void Vector::Normalize(169,4900
+void Vector::Zero(181,5133
+void Vector::One(189,5257
+void Vector::Init(197,5397
+double Vector::Angle(208,5673
+void Vector::RotateVector(217,5940
+Vector& operator+=(operator+=241,6599
+Vector& operator*=(operator*=251,6818
+Vector& operator+(262,7028
+Vector& operator*(275,7281
+bool Vector::Output(287,7510
+ofstream& operator<<(302,7767
+void Vector::Scale(311,7966
+void Vector::Scale(317,8064
+void Vector::Scale(323,8156
+void Vector::Translate(332,8328
+void Vector::MatrixMultiplication(341,8507
+void Vector::InverseMatrixMultiplication(356,8918
+void Vector::LinearCombinationOfVectors(396,10441
+void Vector::Mirror(405,10757
+bool Vector::MakeNormalVector(427,11582
+bool Vector::MakeNormalVector(463,12778
+bool Vector::MakeNormalVector(493,13700
+bool Vector::GetOneNormalVector(512,14285
+double Vector::CutsPlaneAt(558,16093
+bool Vector::LSQdistance(573,16650
+void Vector::AddVector(655,18606
+void Vector::SubtractVector(664,18772
+void Vector::CopyVector(673,18943
+void Vector::AskPosition(684,19258
+bool Vector::SolveSystem(712,20440
+
+
+verbose.cpp,120
+ostream& Verbose::print 7,197
+ostream& operator<<(22,682
+ostream& Binary::print 33,1015
+ostream& operator<<(56,1678
+
+
+graph.cpp,0
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 1ffa216657b5beb85210757a5d292a8b64bd5b79)
+++ src/boundary.cpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -3,7 +3,4 @@
 
 
-// =================spaeter gebrauchte
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol);
-void Find_non_convex_border(Tesselation Tess, molecule mol);
 
 
@@ -29,5 +26,5 @@
 };
 
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line) 
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
   cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
@@ -40,5 +37,5 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryPointSet &a) 
+ostream & operator << (ostream &ost, BoundaryPointSet &a)
 {
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
@@ -73,5 +70,5 @@
 {
   for (int i=0;i<2;i++) {
-    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
     endpoints[i]->lines.erase(Nr);
     LineMap::iterator tester = endpoints[i]->lines.begin();
@@ -85,5 +82,5 @@
 };
 
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
   cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
@@ -92,5 +89,5 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryLineSet &a) 
+ostream & operator << (ostream &ost, BoundaryLineSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
@@ -125,5 +122,5 @@
     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 
+      // and we don't care whether insertion fails
     }
   // set endpoints
@@ -161,5 +158,5 @@
   // get normal vector
   NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x, &endpoints[2]->node->x);
-  
+
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
   if (endpoints[0]->node->x.Projection(&NormalVector) > 0)
@@ -167,5 +164,5 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryTriangleSet &a) 
+ostream & operator << (ostream &ost, BoundaryTriangleSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
@@ -213,5 +210,5 @@
   LineMap LinesOnBoundary;
   TriangleMap TrianglesOnBoundary;
-  
+
   *out << Verbose(1) << "Finding all boundary points." << endl;
   Boundaries *BoundaryPoints = new Boundaries [NDIM];  // first is alpha, second is (r, nr)
@@ -250,5 +247,5 @@
       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) {
@@ -268,8 +265,8 @@
         Walker->x.Output(out);
         *out << endl;
-        double tmp = ProjectedVector.Norm(); 
+        double tmp = ProjectedVector.Norm();
         if (tmp > BoundaryTestPair.first->second.first) {
           BoundaryTestPair.first->second.first = tmp;
-          BoundaryTestPair.first->second.second = Walker; 
+          BoundaryTestPair.first->second.second = Walker;
           *out << Verbose(2) << "Keeping new vector." << endl;
         } else if (tmp == BoundaryTestPair.first->second.first) {
@@ -319,5 +316,5 @@
         }
         // check distance
-        
+
         // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
         {
@@ -328,5 +325,5 @@
   //          SideA.Output(out);
   //          *out << endl;
-          
+
           SideB.CopyVector(&right->second.second->x);
           SideB.ProjectOntoPlane(&AxisVector);
@@ -334,5 +331,5 @@
   //          SideB.Output(out);
   //          *out << endl;
-          
+
           SideC.CopyVector(&left->second.second->x);
           SideC.SubtractVector(&right->second.second->x);
@@ -341,5 +338,5 @@
   //          SideC.Output(out);
   //          *out << endl;
-  
+
           SideH.CopyVector(&runner->second.second->x);
           SideH.ProjectOntoPlane(&AxisVector);
@@ -347,5 +344,5 @@
   //          SideH.Output(out);
   //          *out << endl;
-          
+
           // calculate each length
           double a = SideA.Norm();
@@ -355,9 +352,9 @@
           // calculate the angles
           double alpha = SideA.Angle(&SideH);
-          double beta = SideA.Angle(&SideC); 
+          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(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)) {
@@ -370,5 +367,5 @@
       }
     } while (flag);
-  } 
+  }
   return BoundaryPoints;
 };
@@ -381,5 +378,5 @@
  * \param IsAngstroem whether we have angstroem or atomic units
  * \return NDIM array of the diameters
- */ 
+ */
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
 {
@@ -393,5 +390,5 @@
     *out << Verbose(1) << "Using given boundary points set." << endl;
   }
-  
+
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
@@ -418,5 +415,5 @@
         DistanceVector.SubtractVector(&Neighbour->second.second->x);
         do {  // seek for neighbour pair where it flips
-          OldComponent = DistanceVector.x[Othercomponent]; 
+          OldComponent = DistanceVector.x[Othercomponent];
           Neighbour++;
           if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
@@ -431,5 +428,5 @@
             OtherNeighbour = BoundaryPoints[axis].end();
           OtherNeighbour--;
-          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl; 
+          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
           // now we have found the pair: Neighbour and OtherNeighbour
           OtherVector.CopyVector(&runner->second.second->x);
@@ -466,7 +463,7 @@
  * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
  * \param *mol molecule structure representing the cluster
- * \return determined volume of the cluster in cubed config:GetIsAngstroem() 
+ * \return determined volume of the cluster in cubed config:GetIsAngstroem()
  */
-double VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration, Boundaries *BoundaryPtr, molecule *mol) 
+double VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration, Boundaries *BoundaryPtr, molecule *mol)
 {
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -477,14 +474,14 @@
   double volume = 0.;
   double PyramidVolume = 0.;
-  double G,h; 
+  double G,h;
   Vector x,y;
   double a,b,c;
 
-  Find_non_convex_border(*TesselStruct, *mol);
+  Find_non_convex_border(*TesselStruct, mol);
   /*
   // 1. calculate center of gravity
   *out << endl;
   Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
-  
+
   // 2. translate all points into CoG
   *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
@@ -494,5 +491,5 @@
     Walker->x.Translate(CenterOfGravity);
   }
-  
+
   // 3. Find all points on the boundary
   if (BoundaryPoints == NULL) {
@@ -502,5 +499,5 @@
     *out << Verbose(1) << "Using given boundary points set." << endl;
   }
-  
+
   // 4. fill the boundary point list
   for (int axis=0;axis<NDIM;axis++)
@@ -518,8 +515,8 @@
 //  }
 //  *out << endl;
-  
+
   // 5a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
-  
+
   // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
   TesselStruct->TesselateOnBoundary(out, configuration, mol);
@@ -545,5 +542,5 @@
     PyramidVolume = (1./3.) * G * h;    // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
     *out << Verbose(2) << "Area of triangle is " << G << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " << h << " and the volume is " << PyramidVolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-    volume += PyramidVolume; 
+    volume += PyramidVolume;
   }
   *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
@@ -572,5 +569,5 @@
     for (int i=0;i<mol->AtomCount;i++)
       LookupList[i] = -1;
-    
+
     // print atom coordinates
     *out << Verbose(2) << "The following triangles were created:";
@@ -586,5 +583,5 @@
     for (TriangleMap::iterator runner =  TesselStruct->TrianglesOnBoundary.begin(); runner !=  TesselStruct->TrianglesOnBoundary.end(); runner++) {
       *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 
+      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
     }
     delete[](LookupList);
@@ -595,5 +592,5 @@
   if (BoundaryFreeFlag)
     delete[](BoundaryPoints);
-  
+
   return volume;
 };
@@ -612,5 +609,5 @@
   // transform to PAS
   mol->PrincipalAxisSystem(out, true);
-  
+
   // some preparations beforehand
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -619,5 +616,5 @@
   if (ClusterVolume == 0)
     clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration, BoundaryPoints, mol);
-  else 
+  else
     clustervolume = ClusterVolume;
   double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
@@ -636,7 +633,7 @@
   }
   *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
-  
+
   *out << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass/clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  
+
   // solve cubic polynomial
   *out << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
@@ -647,5 +644,5 @@
     cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
   *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  
+
   double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
   *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
@@ -658,6 +655,6 @@
   } else {
     BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
-    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
-              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2] 
+    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1]
+              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2]
               + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
     BoxLengths.x[2] = minimumvolume - cellvolume;
@@ -669,5 +666,5 @@
       x0 = x2;  // sorted in ascending order
     }
-  
+
     cellvolume = 1;
     for(int i=0;i<NDIM;i++) {
@@ -675,5 +672,5 @@
       cellvolume *= BoxLengths.x[i];
     }
-  
+
     // set new box dimensions
     *out << Verbose(0) << "Translating to box with these boundaries." << endl;
@@ -692,6 +689,6 @@
 Tesselation::Tesselation()
 {
-  PointsOnBoundaryCount = 0; 
-  LinesOnBoundaryCount = 0; 
+  PointsOnBoundaryCount = 0;
+  LinesOnBoundaryCount = 0;
   TrianglesOnBoundaryCount = 0;
 };
@@ -711,5 +708,5 @@
  * \param *out output stream for debugging
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
- */ 
+ */
 void Tesselation::GuessStartingTriangle(ofstream *out)
 {
@@ -722,5 +719,5 @@
   A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
 
-  // with A chosen, take each pair B,C and sort 
+  // with A chosen, take each pair B,C and sort
   if (A != PointsOnBoundary.end()) {
     B = A;
@@ -746,5 +743,5 @@
 //    }
 //    *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
@@ -758,5 +755,5 @@
     PlaneVector.Output(out);
     *out << endl;
-    // 4. loop over all points 
+    // 4. loop over all points
     double sign = 0.;
     PointMap::iterator checker = PointsOnBoundary.begin();
@@ -785,13 +782,13 @@
       tmp = checker->second->node->x.Distance(&A->second->node->x);
       int innerpoint = 0;
-      if ((tmp < A->second->node->x.Distance(&baseline->second.first->second->node->x)) 
+      if ((tmp < A->second->node->x.Distance(&baseline->second.first->second->node->x))
           && (tmp < A->second->node->x.Distance(&baseline->second.second->second->node->x)))
         innerpoint++;
       tmp = checker->second->node->x.Distance(&baseline->second.first->second->node->x);
-      if ((tmp < baseline->second.first->second->node->x.Distance(&A->second->node->x)) 
+      if ((tmp < baseline->second.first->second->node->x.Distance(&A->second->node->x))
           && (tmp < baseline->second.first->second->node->x.Distance(&baseline->second.second->second->node->x)))
         innerpoint++;
       tmp = checker->second->node->x.Distance(&baseline->second.second->second->node->x);
-      if ((tmp < baseline->second.second->second->node->x.Distance(&baseline->second.first->second->node->x)) 
+      if ((tmp < baseline->second.second->second->node->x.Distance(&baseline->second.first->second->node->x))
           && (tmp < baseline->second.second->second->node->x.Distance(&A->second->node->x)))
         innerpoint++;
@@ -815,6 +812,6 @@
     BPS[0] = baseline->second.first->second;
     BPS[1] = A->second;
-    BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
-   
+    BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+
     // 4b3. insert created triangle
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
@@ -825,5 +822,5 @@
       LinesOnBoundaryCount++;
     }
-  
+
     *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
   } else {
@@ -857,7 +854,7 @@
   do {
     flag = false;
-    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) 
+    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
       if (baseline->second->TrianglesCount == 1) {
-        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl; 
+        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
         // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
         SmallestAngle = M_PI;
@@ -888,5 +885,5 @@
         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; 
+        //*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.);
@@ -895,9 +892,9 @@
         *out << endl;
         winner = PointsOnBoundary.end();
-        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) 
+        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++)
           if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
             *out << Verbose(3) << "Target point is " << *(target->second) << ":";
             bool continueflag = true;
-            
+
             VirtualNormalVector.CopyVector(&baseline->second->endpoints[0]->node->x);
             VirtualNormalVector.AddVector(&baseline->second->endpoints[0]->node->x);
@@ -935,5 +932,5 @@
             // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
             continueflag = continueflag && (!(
-                ((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) 
+                ((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
                 && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak))
                ));
@@ -990,10 +987,10 @@
           *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);
-  
+
 };
 
@@ -1021,34 +1018,39 @@
 //====================================================================================================================
 
-
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol)
-{
-  /* d2 ist der Normalenvektor auf dem Dreieck,
-   * d1 ist der Vektor, der normal auf der Kante und d2 steht.
+/*!
+ * This recursive function finds a third point, to form a triangle with two given ones.
+ * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
+ *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
+ *  upon which we operate.
+ *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its id, its general \
+ *  direction and angle into Storage.
+ *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
+ *  with all neighbours of the candidate.
+ */
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule* mol)
+{
+  /* d2 is normal vector on the triangle
+   * d1 is normal on the triangle line, from which we come, as well as on d2.
    */
-  Vector dif_a;
-  Vector dif_b;
-  Vector Chord;
-  Vector AngleCheck;
-  atom *Walker;
-
-  dif_a.CopyVector(&a.x);
-  dif_a.SubtractVector(&Candidate.x);
-  dif_b.CopyVector(&b.x);
-  dif_b.SubtractVector(&Candidate.x);
-  Chord.CopyVector(&a.x);
-  Chord.SubtractVector(&b.x);
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck; // Projection of a difference vector on plane orthogonal on triangle side.
+  atom *Walker; // variable atom point
+
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
   AngleCheck.CopyVector(&dif_a);
-  AngleCheck.ProjectOntoPlane(&Chord);
-
-  //Storage eintrag fuer aktuelles Atom
-  if (Chord.Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  AngleCheck.ProjectOntoPlane(Chord);
+
+  if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
   {
 
     if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
       {
-	Storage[0]=(double)Candidate.nr;
-	Storage[1]=1;
-	Storage[2]=AngleCheck.Angle(d2);
+    	Storage[0]=(double)Candidate->nr;
+    	Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+    	Storage[2]=AngleCheck.Angle(d2);
       }
     else
@@ -1058,5 +1060,5 @@
 	  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
 	  {
-	    Storage[0]=(double)Candidate.nr;
+	    Storage[0]=(double)Candidate->nr;
 	    Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
 	    Storage[2]=AngleCheck.Angle(d2);
@@ -1065,33 +1067,39 @@
   }
 
-  if (n<5)
-    {
-      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
+  if (n<5) // Five is the recursion level threshold.
+    {
+      for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
 	{
-	  while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
-	    {
+		Walker= mol->start; // go through all atoms
+
+	  while (Walker->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr))
+	    { // until atom found which belongs to bond
 	      Walker = Walker->next;
 	    }
 
-	  Find_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
+
+	  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Storage, RADIUS, mol);  //call function again
 	}
     }
 };
 
-
-void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
-{
-  Vector CenterOfLine = Line.endpoints[0]->node->x;
+/*!
+ * this function fins a triangle to a line, adjacent to an existing one.
+ */
+void Tesselation::Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
+{
+	printf("Looking for next suitable triangle \n");
   Vector direction1;
-  Vector direction2;
   Vector helper;
+  Vector Chord;
   atom* Walker;
 
-  double Storage[3];
+  double *Storage;
+  Storage = new double[3];
   Storage[0]=-1.;   // Id must be positive, we see should nothing be done
   Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
   Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
 
-  
+
   helper.CopyVector(&(Line.endpoints[0]->node->x));
   for (int i =0; i<3; i++)
@@ -1114,6 +1122,10 @@
     }
 
-  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &direction1, T.NormalVector, Storage, RADIUS, *mol);
-  
+  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
+  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
+
+printf("Looking for third point candidates for triangle \n");
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, T.NormalVector, Storage, RADIUS, mol);
+
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
   // Next Triangle is Line, atom with number in Storage[0]
@@ -1141,8 +1153,9 @@
   for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
     {
-      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
+
+
+  if ( (LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]))).second);
 	{
-	  LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
-	  LinesOnBoundaryCount++;
+		LinesOnBoundaryCount++;
 	}
     }
@@ -1154,50 +1167,55 @@
       BTS->NormalVector->Scale(-1);
     }
-  
-};
-
-
-void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
-{
+
+};
+
+
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, double Storage[3], molecule* mol)
+{
+	printf("Looking for second point of starting triangle \n");
   int i;
   Vector *AngleCheck;
   atom* Walker;
 
-  AngleCheck->CopyVector(&Candidate.x);
-  AngleCheck->SubtractVector(&a.x);
-  if (AngleCheck->ScalarProduct(&Oben) < Storage[1])
-    {
-      Storage[0]=(double)(Candidate.nr);
-      Storage[1]=AngleCheck->ScalarProduct(&Oben);
+  AngleCheck->CopyVector(&(Candidate->x));
+  AngleCheck->SubtractVector(&(a->x));
+  if (AngleCheck->Angle(&Oben) < Storage[1])
+    {
+	  //printf("Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+      Storage[0]=(double)(Candidate->nr);
+      Storage[1]=AngleCheck->Angle(&Oben);
+      //printf("Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
     };
-
+  printf("%d \n", n);
   if (n<5)
     {
-      for (i = 0; i< mol.NumberOfBondsPerAtom[Candidate.nr]; i++)
+      for (i = 0; i< mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
 	{
-	  Walker = mol.start;
-	  while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
+	  Walker = mol->start;
+	  while (Candidate->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr))
 	    {
 	      Walker = Walker->next;
 	    };
-	  
-	  Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
+
+	  Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Storage, mol);
 	    };
     };
-  
-
-};
-
-
-void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
-{
+
+
+};
+
+
+void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
+{
+	printf("Looking for starting triangle \n");
   int i=0;
-  atom Walker;
-  atom Walker2;
-  atom Walker3;
+  atom* Walker;
+  atom* Walker2;
+  atom* Walker3;
   int max_index[3];
   double max_coordinate[3];
   Vector Oben;
   Vector helper;
+  Vector Chord;
 
   Oben.Zero();
@@ -1209,40 +1227,44 @@
       max_coordinate[i] =-1;
     }
-
-  Walker = *mol.start;
-  while (Walker.next != NULL)
-    {
+printf("Molecule mol is there and has %d Atoms \n", mol->AtomCount);
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
+    	Walker = Walker->next;
       for (i=0; i<3; i++)
-	{
-	  if (Walker.x.x[i]>max_coordinate[i])
+      {
+	  if (Walker->x.x[i] > max_coordinate[i])
 	    {
-	      max_coordinate[i]=Walker.x.x[i];
-	      max_index[i]=Walker.nr;
+	      max_coordinate[i]=Walker->x.x[i];
+	      max_index[i]=Walker->nr;
 	    }
-	}
-    }
-
+      }
+    }
+
+printf("Found starting atom \n");
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-  const int k=0;
-
-  Oben.x[k]=1;
-  Walker = *mol.start;
-  while (Walker.nr != max_index[k])
-    {
-      Walker = *Walker.next;
-    }
-
+  const int k=2;
+
+  Oben.x[k]=1.;
+  Walker = mol->start;
+  Walker = Walker->next;
+  while (Walker->nr != max_index[k])
+    {
+      Walker = Walker->next;
+    }
+printf("%d \n", Walker->nr);
   double Storage[3];
   Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[2]=-10.;  // This will be an angle looking for the third point.
-
-
-  for (i=0; i< mol.NumberOfBondsPerAtom[Walker.nr]; i++)
-    {
-      Walker2 = *mol.start;
-      while (Walker2.nr != (mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom->nr : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr) ) 
+  Storage[1]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[2]=999999.;  // This will be an angle looking for the third point.
+printf("%d \n", mol->NumberOfBondsPerAtom[Walker->nr]);
+
+  for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
+    {
+      Walker2 = mol->start;
+      Walker2 = Walker2->next;
+      while (Walker2->nr != (mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr) )
 	{
-	  Walker2 = *Walker2.next;
+	  Walker2 = Walker2->next;
 	}
 
@@ -1250,38 +1272,45 @@
     }
 
-  Walker2 = *mol.start;
-
-  while (Walker2.nr != int(Storage[0]))
-    {
-      Walker = *Walker.next;
-    }
-  
-  helper.CopyVector(&Walker.x);
-  helper.SubtractVector(&Walker2.x);
+  Walker2 = mol->start;
+
+  while (Walker2->nr != int(Storage[0]))
+    {
+      Walker2 = Walker2->next;
+    }
+
+  helper.CopyVector(&(Walker->x));
+  helper.SubtractVector(&(Walker2->x));
   Oben.ProjectOntoPlane(&helper);
   helper.VectorProduct(&Oben);
-
-       Find_next_suitable_point(Walker, Walker2, *(mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom), 0, &helper, &Oben, Storage, RADIUS, mol);
-  Walker3 = *mol.start;
-  while (Walker3.nr != int(Storage[0]))
-    {
-      Walker3 = *Walker3.next;
+  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
+  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[2]= -10.;  // This will be an angle looking for the third point.
+
+  Chord.CopyVector(&(Walker->x)); // bring into calling function
+  Chord.SubtractVector(&(Walker2->x));
+
+printf("Looking for third point candidates \n");
+       Find_next_suitable_point(Walker, Walker2, (mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][0]->rightatom : mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom), 0, &Chord, &helper, &Oben, Storage, RADIUS, mol);
+  Walker3 = mol->start;
+  while (Walker3->nr != int(Storage[0]))
+    {
+      Walker3 = Walker3->next;
     }
 
   //Starting Triangle is Walker, Walker2, index Storage[0]
 
-  AddPoint(&Walker);
-  AddPoint(&Walker2);
-  AddPoint(&Walker3);
-
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker2);
+  AddPoint(Walker);
+  AddPoint(Walker2);
+  AddPoint(Walker3);
+
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Walker2);
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker3);
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Walker3);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker2);
-  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Walker2);
+  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
@@ -1289,5 +1318,5 @@
   TrianglesOnBoundaryCount++;
 
-  for(int i=0;i<NDIM;i++) 
+  for(int i=0;i<NDIM;i++)
     {
       LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
@@ -1304,14 +1333,14 @@
 
 
-void Find_non_convex_border(Tesselation* Tess, molecule mol) 
-{
+void Find_non_convex_border(Tesselation Tess, molecule* mol)
+{
+	printf("Entering finding of non convex hull. \n");
   const double RADIUS =6;
-  Tess->Find_starting_triangle(mol, RADIUS);
-
-  for (LineMap::iterator baseline = Tess->LinesOnBoundary.begin(); baseline != Tess->LinesOnBoundary.end(); baseline++) 
-    if (baseline->second->TrianglesCount == 1) 
+  Tess.Find_starting_triangle(mol, RADIUS);
+
+  for (LineMap::iterator baseline = Tess.LinesOnBoundary.begin(); baseline != Tess.LinesOnBoundary.end(); baseline++)
+    if (baseline->second->TrianglesCount == 1)
       {
-	Tess->Find_next_suitable_triangle(&mol, baseline->second, baseline->second->triangles.begin()->second, RADIUS); //the line is there, so there is a triangle, but only one.
-
+		Tess.Find_next_suitable_triangle(mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
       }
     else
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 1ffa216657b5beb85210757a5d292a8b64bd5b79)
+++ src/boundary.hpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -82,5 +82,5 @@
     void GuessStartingTriangle(ofstream *out);
     void AddPoint(atom * Walker);
-    void Find_starting_triangle(molecule mol, const double RADIUS);
+    void Find_starting_triangle(molecule* mol, const double RADIUS);
     void Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS);
     
@@ -105,4 +105,6 @@
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
+void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule *mol);
+void Find_non_convex_border(Tesselation Tess, molecule* mol);
 
 
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 1ffa216657b5beb85210757a5d292a8b64bd5b79)
+++ src/molecules.cpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -1,6 +1,6 @@
 /** \file molecules.cpp
- * 
+ *
  * Functions for the class molecule.
- * 
+ *
  */
 
@@ -25,5 +25,5 @@
       sum += (gsl_vector_get(x,j) - (vectors[i])->x[j])*(gsl_vector_get(x,j) - (vectors[i])->x[j]);
   }
-  
+
   return sum;
 };
@@ -34,10 +34,10 @@
  * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
  */
-molecule::molecule(periodentafel *teil) 
-{ 
+molecule::molecule(periodentafel *teil)
+{
   // init atom chain list
-  start = new atom; 
+  start = new atom;
   end = new atom;
-  start->father = NULL; 
+  start->father = NULL;
   end->father = NULL;
   link(start,end);
@@ -46,7 +46,7 @@
   last = new bond(start, end, 1, -1);
   link(first,last);
-  // other stuff 
+  // other stuff
   MDSteps = 0;
-  last_atom = 0; 
+  last_atom = 0;
   elemente = teil;
   AtomCount = 0;
@@ -67,5 +67,5 @@
  * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
  */
-molecule::~molecule() 
+molecule::~molecule()
 {
   if (ListOfBondsPerAtom != NULL)
@@ -78,16 +78,16 @@
   delete(last);
   delete(end);
-  delete(start); 
+  delete(start);
 };
 
 /** Adds given atom \a *pointer from molecule list.
- * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount 
+ * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
  * \param *pointer allocated and set atom
  * \return true - succeeded, false - atom not found in list
  */
 bool molecule::AddAtom(atom *pointer)
-{ 
+{
   if (pointer != NULL) {
-    pointer->sort = &pointer->nr; 
+    pointer->sort = &pointer->nr;
     pointer->nr = last_atom++;  // increase number within molecule
     AtomCount++;
@@ -106,5 +106,5 @@
     return add(pointer, end);
   } else
-    return false; 
+    return false;
 };
 
@@ -115,5 +115,5 @@
  */
 atom * molecule::AddCopyAtom(atom *pointer)
-{ 
+{
   if (pointer != NULL) {
   	atom *walker = new atom();
@@ -122,5 +122,5 @@
     walker->v.CopyVector(&pointer->v); // copy velocity
     walker->FixedIon = pointer->FixedIon;
-    walker->sort = &walker->nr; 
+    walker->sort = &walker->nr;
     walker->nr = last_atom++;  // increase number within molecule
     walker->father = pointer; //->GetTrueFather();
@@ -133,5 +133,5 @@
     return walker;
   } else
-    return NULL; 
+    return NULL;
 };
 
@@ -156,12 +156,12 @@
  *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
  *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
- *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2} 
+ *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
  *    \f]
- *    as the coordination of all three atoms in the coordinate system of these three vectors: 
+ *    as the coordination of all three atoms in the coordinate system of these three vectors:
  *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
- * 
+ *
  * \param *out output stream for debugging
- * \param *Bond pointer to bond between \a *origin and \a *replacement 
- * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin) 
+ * \param *Bond pointer to bond between \a *origin and \a *replacement
+ * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
  * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
  * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
@@ -191,9 +191,9 @@
   InBondvector.SubtractVector(&TopOrigin->x);
   bondlength = InBondvector.Norm();
-   
+
    // is greater than typical bond distance? Then we have to correct periodically
    // the problem is not the H being out of the box, but InBondvector have the wrong direction
-   // due to TopReplacement or Origin being on the wrong side! 
-  if (bondlength > BondDistance) { 
+   // due to TopReplacement or Origin being on the wrong side!
+  if (bondlength > BondDistance) {
 //    *out << Verbose(4) << "InBondvector is: ";
 //    InBondvector.Output(out);
@@ -215,5 +215,5 @@
 //    *out << endl;
   } // periodic correction finished
-  
+
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
@@ -222,5 +222,5 @@
     cerr << Verbose(3) << "ERROR: There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
     return false;
-    BondRescale = bondlength; 
+    BondRescale = bondlength;
   } else {
     if (!IsAngstroem)
@@ -273,5 +273,5 @@
       if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
 //        *out << Verbose(3) << "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane." << endl;
-        
+
         // determine the plane of these two with the *origin
         AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
@@ -286,5 +286,5 @@
       Orthovector1.Normalize();
       //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
-      
+
       // create the two Hydrogens ...
       FirstOtherAtom = new atom();
@@ -318,9 +318,9 @@
         SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle));
       }
-      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance 
+      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
       SecondOtherAtom->x.Scale(&BondRescale);
       //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
       for(int i=NDIM;i--;) { // and make relative to origin atom
-        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i]; 
+        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
         SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
       }
@@ -365,8 +365,8 @@
 //      *out << endl;
       AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
-//      *out << Verbose(3) << "Orthovector2: "; 
+//      *out << Verbose(3) << "Orthovector2: ";
 //      Orthovector2.Output(out);
 //      *out << endl;
-      
+
       // create correct coordination for the three atoms
       alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
@@ -380,10 +380,10 @@
       factors[0] = d;
       factors[1] = f;
-      factors[2] = 0.; 
+      factors[2] = 0.;
       FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
       factors[1] = -0.5*f;
-      factors[2] = g; 
+      factors[2] = g;
       SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
-      factors[2] = -g; 
+      factors[2] = -g;
       ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
 
@@ -437,5 +437,5 @@
  */
 bool molecule::AddXYZFile(string filename)
-{ 
+{
   istringstream *input = NULL;
   int NumberOfAtoms = 0; // atom number in xyz read
@@ -446,5 +446,5 @@
   string line;    // currently parsed line
   double x[3];    // atom coordinates
-  
+
   xyzfile.open(filename.c_str());
   if (!xyzfile)
@@ -454,8 +454,8 @@
   input = new istringstream(line);
   *input >> NumberOfAtoms;
-  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl; 
+  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
   getline(xyzfile,line,'\n'); // Read comment
   cout << Verbose(1) << "Comment: " << line << endl;
-  
+
   if (MDSteps == 0) // no atoms yet present
     MDSteps++;
@@ -491,5 +491,5 @@
   xyzfile.close();
   delete(input);
-  return true; 
+  return true;
 };
 
@@ -503,5 +503,5 @@
   atom *LeftAtom = NULL, *RightAtom = NULL;
   atom *Walker = NULL;
-  
+
   // copy all atoms
   Walker = start;
@@ -510,5 +510,5 @@
     CurrentAtom = copy->AddCopyAtom(Walker);
   }
-  
+
   // copy all bonds
   bond *Binder = first;
@@ -534,5 +534,5 @@
       copy->NoCyclicBonds++;
     NewBond->Type = Binder->Type;
-  }  
+  }
   // correct fathers
   Walker = copy->start;
@@ -551,5 +551,5 @@
     copy->CreateListOfBondsPerAtom((ofstream *)&cout);
   }
-  
+
   return copy;
 };
@@ -576,5 +576,5 @@
 
 /** Remove bond from bond chain list.
- * \todo Function not implemented yet 
+ * \todo Function not implemented yet
  * \param *pointer bond pointer
  * \return true - bound found and removed, false - bond not found/removed
@@ -588,5 +588,5 @@
 
 /** Remove every bond from bond chain list that atom \a *BondPartner is a constituent of.
- * \todo Function not implemented yet 
+ * \todo Function not implemented yet
  * \param *BondPartner atom to be removed
  * \return true - bounds found and removed, false - bonds not found/removed
@@ -621,5 +621,5 @@
   Vector *min = new Vector;
   Vector *max = new Vector;
-  
+
   // gather min and max for each axis
   ptr = start->next;  // start at first in list
@@ -667,5 +667,5 @@
 {
   Vector *min = new Vector;
-  
+
 //  *out << Verbose(3) << "Begin of CenterEdge." << endl;
   atom *ptr = start->next;  // start at first in list
@@ -683,5 +683,5 @@
       }
     }
-//    *out << Verbose(4) << "Maximum is "; 
+//    *out << Verbose(4) << "Maximum is ";
 //    max->Output(out);
 //    *out << ", Minimum is ";
@@ -691,8 +691,8 @@
     max->AddVector(min);
     Translate(min);
-  } 
+  }
   delete(min);
 //  *out << Verbose(3) << "End of CenterEdge." << endl;
-}; 
+};
 
 /** Centers the center of the atoms at (0,0,0).
@@ -704,18 +704,18 @@
   int Num = 0;
   atom *ptr = start->next;  // start at first in list
- 
+
   for(int i=NDIM;i--;) // zero center vector
     center->x[i] = 0.;
-    
+
   if (ptr != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
       ptr = ptr->next;
       Num++;
-      center->AddVector(&ptr->x);       
+      center->AddVector(&ptr->x);
     }
     center->Scale(-1./Num); // divide through total number (and sign for direction)
     Translate(center);
   }
-}; 
+};
 
 /** Returns vector pointing to center of gravity.
@@ -729,5 +729,5 @@
   Vector tmp;
   double Num = 0;
-  
+
   a->Zero();
 
@@ -737,5 +737,5 @@
       Num += 1.;
       tmp.CopyVector(&ptr->x);
-      a->AddVector(&tmp);       
+      a->AddVector(&tmp);
     }
     a->Scale(-1./Num); // divide through total mass (and sign for direction)
@@ -757,5 +757,5 @@
 	Vector tmp;
   double Num = 0;
-	
+
 	a->Zero();
 
@@ -766,5 +766,5 @@
       tmp.CopyVector(&ptr->x);
       tmp.Scale(ptr->type->mass);  // scale by mass
-      a->AddVector(&tmp);       
+      a->AddVector(&tmp);
     }
     a->Scale(-1./Num); // divide through total mass (and sign for direction)
@@ -789,5 +789,5 @@
     Translate(center);
   }
-}; 
+};
 
 /** Scales all atoms by \a *factor.
@@ -803,8 +803,8 @@
       Trajectories[ptr].R.at(j).Scale(factor);
     ptr->x.Scale(factor);
-  }      
-};
-
-/** Translate all atoms by given vector. 
+  }
+};
+
+/** Translate all atoms by given vector.
  * \param trans[] translation vector.
  */
@@ -818,8 +818,8 @@
       Trajectories[ptr].R.at(j).Translate(trans);
     ptr->x.Translate(trans);
-  }      
-};
-
-/** Mirrors all atoms against a given plane. 
+  }
+};
+
+/** Mirrors all atoms against a given plane.
  * \param n[] normal vector of mirror plane.
  */
@@ -833,5 +833,5 @@
       Trajectories[ptr].R.at(j).Mirror(n);
     ptr->x.Mirror(n);
-  }      
+  }
 };
 
@@ -847,5 +847,5 @@
   bool flag;
   Vector Testvector, Translationvector;
-  
+
   do {
     Center.Zero();
@@ -863,5 +863,5 @@
           if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j]; 
+              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
@@ -879,5 +879,5 @@
         cout << Verbose(1) << "vector is: ";
         Testvector.Output((ofstream *)&cout);
-        cout << endl;     
+        cout << endl;
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
@@ -892,5 +892,5 @@
             cout << Verbose(1) << "Hydrogen vector is: ";
             Testvector.Output((ofstream *)&cout);
-            cout << endl;     
+            cout << endl;
           }
         }
@@ -914,9 +914,9 @@
 
 	CenterGravity(out, CenterOfGravity);
-	
-	// reset inertia tensor	
+
+	// reset inertia tensor
 	for(int i=0;i<NDIM*NDIM;i++)
 		InertiaTensor[i] = 0.;
-	
+
 	// sum up inertia tensor
 	while (ptr->next != end) {
@@ -943,5 +943,5 @@
 	}
 	*out << endl;
-	
+
 	// diagonalize to determine principal axis system
 	gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
@@ -952,13 +952,13 @@
 	gsl_eigen_symmv_free(T);
 	gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
-	
+
 	for(int i=0;i<NDIM;i++) {
 		*out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
 		*out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
 	}
-	
+
 	// check whether we rotate or not
 	if (DoRotate) {
-	  *out << Verbose(1) << "Transforming molecule into PAS ... "; 
+	  *out << Verbose(1) << "Transforming molecule into PAS ... ";
 	  // the eigenvectors specify the transformation matrix
 	  ptr = start;
@@ -972,8 +972,8 @@
 
 	  // summing anew for debugging (resulting matrix has to be diagonal!)
-	  // reset inertia tensor 
+	  // reset inertia tensor
     for(int i=0;i<NDIM*NDIM;i++)
       InertiaTensor[i] = 0.;
-    
+
     // sum up inertia tensor
     ptr = start;
@@ -1002,5 +1002,5 @@
     *out << endl;
 	}
-	
+
 	// free everything
 	delete(CenterOfGravity);
@@ -1031,5 +1031,5 @@
 
   CountElements();  // make sure ElementsInMolecule is up to date
-  
+
   // check file
   if (input == NULL) {
@@ -1087,5 +1087,5 @@
               Trajectories[walker].U.at(MDSteps).x[d] += 0.5*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d]+Trajectories[walker].F.at(MDSteps).x[d])/IonMass;
             }
-//            cout << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
+//            cout << "Integrated position&velocity of step " << (MDSteps) << ": (";
 //            for (int d=0;d<NDIM;d++)
 //              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
@@ -1121,10 +1121,10 @@
   MDSteps++;
 
-  
+
   // exit
   return true;
 };
 
-/** Align all atoms in such a manner that given vector \a *n is along z axis. 
+/** Align all atoms in such a manner that given vector \a *n is along z axis.
  * \param n[] alignment vector.
  */
@@ -1145,12 +1145,12 @@
     ptr = ptr->next;
     tmp = ptr->x.x[0];
-    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2]; 
+    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     for (int j=0;j<MDSteps;j++) {
       tmp = Trajectories[ptr].R.at(j).x[0];
-      Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2]; 
+      Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
     }
-  }      
+  }
   // rotate n vector
   tmp = n->x[0];
@@ -1160,5 +1160,5 @@
   n->Output((ofstream *)&cout);
   cout << endl;
-  
+
   // rotate on z-y plane
   ptr = start;
@@ -1168,17 +1168,17 @@
     ptr = ptr->next;
     tmp = ptr->x.x[1];
-    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2]; 
+    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     for (int j=0;j<MDSteps;j++) {
       tmp = Trajectories[ptr].R.at(j).x[1];
-      Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2]; 
+      Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
     }
-  }      
+  }
   // rotate n vector (for consistency check)
   tmp = n->x[1];
   n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
   n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
-  
+
   cout << Verbose(1) << "alignment vector after second rotation: ";
   n->Output((ofstream *)&cout);
@@ -1191,10 +1191,10 @@
  * \return true - succeeded, false - atom not found in list
  */
-bool molecule::RemoveAtom(atom *pointer)  
-{ 
+bool molecule::RemoveAtom(atom *pointer)
+{
   if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
     ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
   else
-    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl; 
+    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
@@ -1206,7 +1206,7 @@
  * \return true - succeeded, false - atom not found in list
  */
-bool molecule::CleanupMolecule() 
-{ 
-  return (cleanup(start,end) && cleanup(first,last)); 
+bool molecule::CleanupMolecule()
+{
+  return (cleanup(start,end) && cleanup(first,last));
 };
 
@@ -1222,5 +1222,5 @@
   } else {
     cout << Verbose(0) << "Atom not found in list." << endl;
-    return NULL;  
+    return NULL;
   }
 };
@@ -1271,5 +1271,5 @@
   struct lsq_params *par = (struct lsq_params *)params;
   atom *ptr = par->mol->start;
-  
+
   // initialize vectors
   a.x[0] = gsl_vector_get(x,0);
@@ -1301,5 +1301,5 @@
 {
     int np = 6;
-    
+
    const gsl_multimin_fminimizer_type *T =
      gsl_multimin_fminimizer_nmsimplex;
@@ -1307,19 +1307,19 @@
    gsl_vector *ss;
    gsl_multimin_function minex_func;
- 
+
    size_t iter = 0, i;
    int status;
    double size;
- 
+
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
- 
+
    /* Set all step sizes to 1 */
    gsl_vector_set_all (ss, 1.0);
- 
+
    /* Starting point */
    par->x = gsl_vector_alloc (np);
    par->mol = this;
- 
+
    gsl_vector_set (par->x, 0, 0.0);  // offset
    gsl_vector_set (par->x, 1, 0.0);
@@ -1328,29 +1328,29 @@
    gsl_vector_set (par->x, 4, 0.0);
    gsl_vector_set (par->x, 5, 1.0);
- 
+
    /* Initialize method and iterate */
    minex_func.f = &LeastSquareDistance;
    minex_func.n = np;
    minex_func.params = (void *)par;
- 
+
    s = gsl_multimin_fminimizer_alloc (T, np);
    gsl_multimin_fminimizer_set (s, &minex_func, par->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)
          {
            printf ("converged to minimum at\n");
          }
- 
+
        printf ("%5d ", (int)iter);
        for (i = 0; i < (size_t)np; i++)
@@ -1361,5 +1361,5 @@
      }
    while (status == GSL_CONTINUE && iter < 100);
- 
+
   for (i=0;i<(size_t)np;i++)
     gsl_vector_set(par->x, i, gsl_vector_get(s->x, i));
@@ -1378,5 +1378,5 @@
   int ElementNo, AtomNo;
   CountElements();
-  
+
   if (out == NULL) {
     return false;
@@ -1413,5 +1413,5 @@
   int ElementNo, AtomNo;
   CountElements();
-  
+
   if (out == NULL) {
     return false;
@@ -1460,5 +1460,5 @@
   atom *Walker = start;
   while (Walker->next != end) {
-    Walker = Walker->next; 
+    Walker = Walker->next;
 #ifdef ADDHYDROGEN
     if (Walker->type->Z != 1) {   // regard only non-hydrogen
@@ -1491,5 +1491,5 @@
   int No = 0;
   time_t now;
-  
+
   now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
   walker = start;
@@ -1520,5 +1520,5 @@
   int No = 0;
   time_t now;
-  
+
   now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
   walker = start;
@@ -1563,5 +1563,5 @@
 	      Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
 	      if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
-	        NoNonHydrogen++; 
+	        NoNonHydrogen++;
 	      Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
 	      Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
@@ -1571,5 +1571,5 @@
 	    }
     } else
-    	*out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl; 
+    	*out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
   }
 };
@@ -1583,5 +1583,5 @@
   	ElementsInMolecule[i] = 0;
 	ElementCount = 0;
-	
+
   atom *walker = start;
   while (walker->next != end) {
@@ -1619,5 +1619,5 @@
     Binder = Binder->next;
     if (Binder->Cyclic)
-      No++;    
+      No++;
   }
   delete(BackEdgeStack);
@@ -1679,5 +1679,5 @@
  * Generally, we use the CSD approach to bond recognition, that is the the distance
  * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
- * a threshold t = 0.4 Angstroem. 
+ * a threshold t = 0.4 Angstroem.
  * To make it O(N log N) the function uses the linked-cell technique as follows:
  * The procedure is step-wise:
@@ -1694,6 +1694,50 @@
  * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
  */
+void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)// ofstream* out, double bonddistance, bool IsAngstroem)
+{
+
+	// 1 We will parse bonds out of the dbond file created by tremolo.
+			int atom1, atom2, temp;
+			atom *Walker, *OtherWalker;
+
+	        if (!input)
+	        {
+	        	cout << Verbose(1) << "Opening silica failed \n";
+	        };
+
+			*input >> ws >> atom1;
+			*input >> ws >> atom2;
+	        cout << Verbose(1) << "Scanning file\n";
+	        while (!input->eof()) // Check whether we read 2 items
+	        {
+				*input >> ws >> atom1;
+				*input >> ws >> atom2;
+	        	if(atom2<atom1) //Sort indizes of atoms in order
+	        	{
+	        		temp=atom1;
+	        		atom1=atom2;
+	        		atom2=temp;
+	        	};
+
+	        	Walker=start;
+	        	while(Walker-> nr != atom1) // Find atom corresponding to first index
+	        	{
+	        		Walker = Walker->next;
+	        	};
+	        	OtherWalker = Walker->next;
+	        	while(OtherWalker->nr != atom2) // Find atom corresponding to second index
+	        	{
+	        		OtherWalker= OtherWalker->next;
+	        	};
+	        	AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
+	        	BondCount++; //Increase bond count of molecule
+	        }
+
+	        CreateListOfBondsPerAtom(out);
+
+};
 void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
 {
+
   atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
   int No, NoBonds, CandidateBondNo;
@@ -1704,5 +1748,5 @@
   Vector x;
   int FalseBondDegree = 0;
- 
+
   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
   *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
@@ -1711,5 +1755,5 @@
     cleanup(first,last);
   }
-	
+
   // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
   CountAtoms(out);
@@ -1730,5 +1774,5 @@
     for (int i=NumberCells;i--;)
       CellList[i] = NULL;
-  
+
     // 2b. put all atoms into its corresponding list
     Walker = start;
@@ -1751,10 +1795,11 @@
       if (CellList[index] == NULL)  // allocate molecule if not done
         CellList[index] = new molecule(elemente);
-      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference 
-      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl; 
+      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
+      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
     }
     //for (int i=0;i<NumberCells;i++)
       //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
-      
+
+
     // 3a. go through every cell
     for (N[0]=divisor[0];N[0]--;)
@@ -1769,5 +1814,5 @@
               Walker = Walker->next;
               //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
-              // 3c. check for possible bond between each atom in this and every one in the 27 cells 
+              // 3c. check for possible bond between each atom in this and every one in the 27 cells
               for (n[0]=-1;n[0]<=1;n[0]++)
                 for (n[1]=-1;n[1]<=1;n[1]++)
@@ -1801,4 +1846,7 @@
           }
         }
+
+
+
     // 4. free the cell again
     for (int i=NumberCells;i--;)
@@ -1807,8 +1855,8 @@
       }
     Free((void **)&CellList, "molecule::CreateAdjacencyList - ** CellList");
-    
+
     // create the adjacency list per atom
     CreateListOfBondsPerAtom(out);
-  		
+
     // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
     // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
@@ -1859,5 +1907,5 @@
 			*out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
 	  *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
-		
+
 	  // output bonds for debugging (if bond chain list was correctly installed)
 	  *out << Verbose(1) << endl << "From contents of bond chain list:";
@@ -1869,8 +1917,11 @@
     *out << endl;
   } else
-  	*out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl; 
+  	*out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
   Free((void **)&matrix, "molecule::CreateAdjacencyList: *matrix");
-};
+
+};
+
+
 
 /** Performs a Depth-First search on this molecule.
@@ -1893,7 +1944,7 @@
   bond *Binder = NULL;
   bool BackStepping = false;
-  
+
   *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
-  
+
   ResetAllBondsToUnused();
   ResetAllAtomNumbers();
@@ -1908,5 +1959,5 @@
     LeafWalker->Leaf = new molecule(elemente);
     LeafWalker->Leaf->AddCopyAtom(Root);
-    
+
     OldGraphNr = CurrentGraphNr;
     Walker = Root;
@@ -1919,5 +1970,5 @@
           AtomStack->Push(Walker);
           CurrentGraphNr++;
-        }      
+        }
         do { // (3) if Walker has no unused egdes, go to (5)
           BackStepping = false; // reset backstepping flag for (8)
@@ -1953,10 +2004,10 @@
           Binder = NULL;
       } while (1);  // (2)
-      
+
       // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
       if ((Walker == Root) && (Binder == NULL))
         break;
-        
-      // (5) if Ancestor of Walker is ... 
+
+      // (5) if Ancestor of Walker is ...
       *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
       if (Walker->Ancestor->GraphNr != Root->GraphNr) {
@@ -2001,5 +2052,5 @@
         } while (OtherAtom != Walker);
         ComponentNumber++;
-    
+
         // (11) Root is separation vertex,  set Walker to Root and go to (4)
         Walker = Root;
@@ -2014,5 +2065,5 @@
 
     // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
-    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;    
+    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
     LeafWalker->Leaf->Output(out);
     *out << endl;
@@ -2022,5 +2073,5 @@
       //*out << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
       if (Root->GraphNr != -1) // if already discovered, step on
-        Root = Root->next; 
+        Root = Root->next;
     }
   }
@@ -2044,5 +2095,5 @@
     *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
   }
-  
+
   *out << Verbose(1) << "Final graph info for each bond is:" << endl;
   Binder = first;
@@ -2055,5 +2106,5 @@
     *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
     OutputComponentNumber(out, Binder->rightatom);
-    *out << ">." << endl; 
+    *out << ">." << endl;
     if (Binder->Cyclic) // cyclic ??
       *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
@@ -2070,5 +2121,5 @@
  * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
  * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
- * as cyclic and print out the cycles. 
+ * as cyclic and print out the cycles.
  * \param *out output stream for debugging
  * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
@@ -2081,5 +2132,5 @@
   int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
   enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
-  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring 
+  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring
   class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
@@ -2093,8 +2144,8 @@
     ColorList[i] = white;
   }
-  
+
   *out << Verbose(1) << "Back edge list - ";
   BackEdgeStack->Output(out);
-  
+
   *out << Verbose(1) << "Analysing cycles ... " << endl;
   NumCycles = 0;
@@ -2102,7 +2153,7 @@
     BackEdge = BackEdgeStack->PopFirst();
     // this is the target
-    Root = BackEdge->leftatom;  
+    Root = BackEdge->leftatom;
     // this is the source point
-    Walker = BackEdge->rightatom; 
+    Walker = BackEdge->rightatom;
     ShortestPathList[Walker->nr] = 0;
     BFSStack->ClearStack();  // start with empty BFS stack
@@ -2118,5 +2169,5 @@
         if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
           OtherAtom = Binder->GetOtherAtom(Walker);
-#ifdef ADDHYDROGEN          
+#ifdef ADDHYDROGEN
           if (OtherAtom->type->Z != 1) {
 #endif
@@ -2127,5 +2178,5 @@
               PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
               ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-              *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+              *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
               //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
                 *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
@@ -2137,5 +2188,5 @@
             if (OtherAtom == Root)
               break;
-#ifdef ADDHYDROGEN          
+#ifdef ADDHYDROGEN
           } else {
             *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
@@ -2175,5 +2226,5 @@
       }
     } while ((!BFSStack->IsEmpty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
-    
+
     if (OtherAtom == Root) {
       // now climb back the predecessor list and thus find the cycle members
@@ -2203,5 +2254,5 @@
       *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
     }
-    
+
     // now clean the lists
     while (!TouchedStack->IsEmpty()){
@@ -2213,9 +2264,9 @@
   }
   if (MinRingSize != -1) {
-    // go over all atoms 
+    // go over all atoms
     Root = start;
     while(Root->next != end) {
       Root = Root->next;
-      
+
       if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
         Walker = Root;
@@ -2254,7 +2305,7 @@
           }
           ColorList[Walker->nr] = black;
-          //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
+          //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
         }
-    
+
         // now clean the lists
         while (!TouchedStack->IsEmpty()){
@@ -2305,5 +2356,5 @@
 void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
 {
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++) 
+  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
     *out << vertex->ComponentNr[i] << "  ";
 };
@@ -2383,5 +2434,5 @@
 {
   int c = 0;
-  int FragmentCount; 
+  int FragmentCount;
   // get maximum bond degree
   atom *Walker = start;
@@ -2393,10 +2444,10 @@
   *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
   return FragmentCount;
-};  
+};
 
 /** Scans a single line for number and puts them into \a KeySet.
  * \param *out output stream for debugging
  * \param *buffer buffer to scan
- * \param &CurrentSet filled KeySet on return 
+ * \param &CurrentSet filled KeySet on return
  * \return true - at least one valid atom id parsed, false - CurrentSet is empty
  */
@@ -2406,5 +2457,5 @@
   int AtomNr;
   int status = 0;
-  
+
   line.str(buffer);
   while (!line.eof()) {
@@ -2442,9 +2493,9 @@
   double TEFactor;
   char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
-  
+
   if (FragmentList == NULL) { // check list pointer
     FragmentList = new Graph;
   }
-  
+
   // 1st pass: open file and read
   *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
@@ -2475,5 +2526,5 @@
     status = false;
   }
-  
+
   // 2nd pass: open TEFactors file and read
   *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
@@ -2487,5 +2538,5 @@
         InputFile >> TEFactor;
         (*runner).second.second = TEFactor;
-        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl; 
+        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
       } else {
         status = false;
@@ -2528,5 +2579,5 @@
   if(output != NULL) {
     for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
-      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) { 
+      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
         if (sprinter != (*runner).first.begin())
           output << "\t";
@@ -2594,5 +2645,5 @@
     status = false;
   }
-  
+
   return status;
 };
@@ -2603,5 +2654,5 @@
  * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
  * \return true - structure is equal, false - not equivalence
- */ 
+ */
 bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
 {
@@ -2610,5 +2661,5 @@
   bool status = true;
   char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
-  
+
   filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
   File.open(filename.str().c_str(), ios::out);
@@ -2669,5 +2720,5 @@
   *out << endl;
   Free((void **)&buffer, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
-  
+
   return status;
 };
@@ -2691,5 +2742,5 @@
   for(int i=AtomCount;i--;)
     AtomMask[i] = false;
-  
+
   if (Order < 0) { // adaptive increase of BondOrder per site
     if (AtomMask[AtomCount] == true)  // break after one step
@@ -2731,5 +2782,5 @@
           line >> ws >> Value; // skip time entry
           line >> ws >> Value;
-          No -= 1;  // indices start at 1 in file, not 0 
+          No -= 1;  // indices start at 1 in file, not 0
           //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
 
@@ -2740,8 +2791,8 @@
             // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
             pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
-            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first; 
+            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
             if (!InsertedElement.second) { // this root is already present
-              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term 
-                //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase) 
+              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
+                //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
                 {  // if value is smaller, update value and order
                 (*PresentItem).second.first = fabs(Value);
@@ -2781,9 +2832,9 @@
         Walker = FindAtom(No);
         //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
-          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl; 
+          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
           AtomMask[No] = true;
           status = true;
         //} else
-          //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl; 
+          //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
       }
       // close and done
@@ -2819,5 +2870,5 @@
     if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
       status = true;
-      
+
     if (!status) {
       if (Order == 0)
@@ -2827,5 +2878,5 @@
     }
   }
-  
+
   // print atom mask for debugging
   *out << "              ";
@@ -2836,5 +2887,5 @@
     *out << (AtomMask[i] ? "t" : "f");
   *out << endl;
- 
+
   return status;
 };
@@ -2850,5 +2901,5 @@
   int AtomNo = 0;
   atom *Walker = NULL;
-  
+
   if (SortIndex != NULL) {
     *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
@@ -2908,6 +2959,6 @@
   atom **ListOfAtoms = NULL;
   atom ***ListOfLocalAtoms = NULL;
-  bool *AtomMask = NULL; 
-  
+  bool *AtomMask = NULL;
+
   *out << endl;
 #ifdef ADDHYDROGEN
@@ -2918,7 +2969,7 @@
 
   // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
-  
+
   // ===== 1. Check whether bond structure is same as stored in files ====
-  
+
   // fill the adjacency list
   CreateListOfBondsPerAtom(out);
@@ -2926,7 +2977,7 @@
   // create lookup table for Atom::nr
   FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);
-  
+
   // === compare it with adjacency file ===
-  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms); 
+  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
   Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
 
@@ -2957,5 +3008,5 @@
     delete(LocalBackEdgeStack);
   }
-  
+
   // ===== 3. if structure still valid, parse key set file and others =====
   FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
@@ -2963,7 +3014,7 @@
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
   FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
-  
-  // =================================== Begin of FRAGMENTATION =============================== 
-  // ===== 6a. assign each keyset to its respective subgraph ===== 
+
+  // =================================== Begin of FRAGMENTATION ===============================
+  // ===== 6a. assign each keyset to its respective subgraph =====
   Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
 
@@ -2976,5 +3027,5 @@
     FragmentationToDo = FragmentationToDo || CheckOrder;
     AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
-    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) ===== 
+    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
     Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
 
@@ -3000,5 +3051,5 @@
   delete(ParsedFragmentList);
   delete[](MinimumRingSize);
-  
+
 
   // ==================================== End of FRAGMENTATION ============================================
@@ -3006,5 +3057,5 @@
   // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
   Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
-  
+
   // free subgraph memory again
   FragmentCounter = 0;
@@ -3031,10 +3082,10 @@
     }
     *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
-    
+
     // ===== 9. Save fragments' configuration and keyset files et al to disk ===
     if (BondFragments->NumberOfMolecules != 0) {
       // create the SortIndex from BFS labels to order in the config file
       CreateMappingLabelsToConfigSequence(out, SortIndex);
-      
+
       *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
       if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
@@ -3042,30 +3093,30 @@
       else
         *out << Verbose(1) << "Some config writing failed." << endl;
-  
+
       // store force index reference file
       BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
-      
-      // store keysets file 
+
+      // store keysets file
       StoreKeySetFile(out, TotalGraph, configuration->configpath);
-  
-      // store Adjacency file 
+
+      // store Adjacency file
       StoreAdjacencyToFile(out, configuration->configpath);
-  
+
       // store Hydrogen saturation correction file
       BondFragments->AddHydrogenCorrection(out, configuration->configpath);
-      
+
       // store adaptive orders into file
       StoreOrderAtSiteFile(out, configuration->configpath);
-      
+
       // restore orbital and Stop values
       CalculateOrbitals(*configuration);
-      
+
       // free memory for bond part
       *out << Verbose(1) << "Freeing bond memory" << endl;
       delete(FragmentList); // remove bond molecule from memory
-      Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex"); 
+      Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
     } else
       *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
-  //} else 
+  //} else
   //  *out << Verbose(1) << "No fragments to store." << endl;
   *out << Verbose(0) << "End of bond fragmentation." << endl;
@@ -3093,8 +3144,8 @@
   atom *Walker = NULL, *OtherAtom = NULL;
   ReferenceStack->Push(Binder);
-  
+
   do {  // go through all bonds and push local ones
     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
-    if (Walker != NULL) // if this Walker exists in the subgraph ... 
+    if (Walker != NULL) // if this Walker exists in the subgraph ...
     	for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
       	OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
@@ -3109,5 +3160,5 @@
     ReferenceStack->Push(Binder);
   } while (FirstBond != Binder);
-  
+
   return status;
 };
@@ -3185,5 +3236,5 @@
       Walker->AdaptiveOrder = OrderArray[Walker->nr];
       Walker->MaxOrder = MaxArray[Walker->nr];
-      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl; 
+      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl;
     }
     file.close();
@@ -3196,5 +3247,5 @@
   Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
   Free((void **)&MaxArray, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
-  
+
   *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
   return status;
@@ -3254,5 +3305,5 @@
   Walker = start;
   while (Walker->next != end) {
-    Walker = Walker->next; 
+    Walker = Walker->next;
     *out << Verbose(4) << "Atom " << Walker->Name << "/" << Walker->nr << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
     TotalDegree = 0;
@@ -3262,5 +3313,5 @@
     }
     *out << " -- TotalDegree: " << TotalDegree << endl;
-  }      
+  }
   *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
 };
@@ -3268,5 +3319,5 @@
 /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
  * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
- * white and putting into queue. 
+ * white and putting into queue.
  * \param *out output stream for debugging
  * \param *Mol Molecule class to add atoms to
@@ -3277,5 +3328,5 @@
  * \param BondOrder maximum distance for vertices to add
  * \param IsAngstroem lengths are in angstroem or bohrradii
- */  
+ */
 void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
 {
@@ -3303,5 +3354,5 @@
   }
   ShortestPathList[Root->nr] = 0;
-  
+
   // and go on ... Queue always contains all lightgray vertices
   while (!AtomStack->IsEmpty()) {
@@ -3311,5 +3362,5 @@
     // followed by n+1 till top of stack.
     Walker = AtomStack->PopFirst(); // pop oldest added
-    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl; 
+    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
     for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
       Binder = ListOfBondsPerAtom[Walker->nr][i];
@@ -3318,9 +3369,9 @@
         *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
         if (ColorList[OtherAtom->nr] == white) {
-          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem) 
+          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
             ColorList[OtherAtom->nr] = lightgray;
           PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
           ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-          *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+          *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
           if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
             *out << Verbose(3);
@@ -3370,5 +3421,5 @@
           // This has to be a cyclic bond, check whether it's present ...
           if (AddedBondList[Binder->nr] == NULL) {
-            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) { 
+            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
               AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
               AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
@@ -3385,5 +3436,5 @@
     }
     ColorList[Walker->nr] = black;
-    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
+    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
   }
   Free((void **)&PredecessorList, "molecule::BreadthFirstSearchAdd: **PredecessorList");
@@ -3409,10 +3460,10 @@
 
   *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
-  
+
   // reset parent list
   *out << Verbose(3) << "Resetting ParentList." << endl;
   for (int i=Father->AtomCount;i--;)
     ParentList[i] = NULL;
-  
+
   // fill parent list with sons
   *out << Verbose(3) << "Filling Parent List." << endl;
@@ -3455,5 +3506,5 @@
  * \param *&Leaf KeySet to look through
  * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
- * \param index of the atom suggested for removal 
+ * \param index of the atom suggested for removal
  */
 int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
@@ -3461,5 +3512,5 @@
   atom *Runner = NULL;
   int SP, Removal;
-  
+
   *out << Verbose(2) << "Looking for removal candidate." << endl;
   SP = -1; //0;  // not -1, so that Root is never removed
@@ -3479,7 +3530,7 @@
 /** Stores a fragment from \a KeySet into \a molecule.
  * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
- * molecule and adds missing hydrogen where bonds were cut. 
+ * molecule and adds missing hydrogen where bonds were cut.
  * \param *out output stream for debugging messages
- * \param &Leaflet pointer to KeySet structure 
+ * \param &Leaflet pointer to KeySet structure
  * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return pointer to constructed molecule
@@ -3492,10 +3543,10 @@
   bool LonelyFlag = false;
   int size;
-  
+
 //  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
-  
+
   Leaf->BondDistance = BondDistance;
   for(int i=NDIM*2;i--;)
-    Leaf->cell_size[i] = cell_size[i]; 
+    Leaf->cell_size[i] = cell_size[i];
 
   // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
@@ -3510,5 +3561,5 @@
     size++;
   }
-  
+
   // create the bonds between all: Make it an induced subgraph and add hydrogen
 //  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
@@ -3520,21 +3571,21 @@
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
-      for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father 
+      for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
         OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
 //        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
         if (SonList[OtherFather->nr] != NULL) {
-//          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl; 
+//          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
           if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
 //            *out << Verbose(3) << "Adding Bond: ";
-//            *out << 
+//            *out <<
             Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
 //            *out << "." << endl;
             //NumBonds[Runner->nr]++;
-          } else { 
+          } else {
 //            *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
           }
           LonelyFlag = false;
         } else {
-//          *out << ", who has no son in this fragment molecule." << endl; 
+//          *out << ", who has no son in this fragment molecule." << endl;
 #ifdef ADDHYDROGEN
           //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
@@ -3554,5 +3605,5 @@
     while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
       Runner = Runner->next;
-#endif       
+#endif
   }
   Leaf->CreateListOfBondsPerAtom(out);
@@ -3587,5 +3638,5 @@
   StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
   StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
-  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL; 
+  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
   MoleculeListClass *FragmentList = NULL;
   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
@@ -3637,5 +3688,5 @@
 		// clear snake stack
 	  SnakeStack->ClearStack();
-    //SnakeStack->TestImplementation(out, start->next);    
+    //SnakeStack->TestImplementation(out, start->next);
 
     ///// INNER LOOP ////////////
@@ -3658,5 +3709,5 @@
         }
         if (ShortestPathList[Walker->nr] == -1) {
-          ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1; 
+          ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1;
           TouchedStack->Push(Walker); // mark every atom for lists cleanup later, whose shortest path has been changed
         }
@@ -3696,10 +3747,10 @@
     				OtherAtom = Binder->GetOtherAtom(Walker);
             if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
-              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl; 
+              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
               //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
             } else { // otherwise check its colour and element
       				if (
 #ifdef ADDHYDROGEN
-              (OtherAtom->type->Z != 1) && 
+              (OtherAtom->type->Z != 1) &&
 #endif
                     (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
@@ -3711,11 +3762,11 @@
                 //} else {
                 //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
-                //} 
+                //}
                 Walker = OtherAtom;
                 break;
               } else {
-                if (OtherAtom->type->Z == 1) 
+                if (OtherAtom->type->Z == 1)
                   *out << "Links to a hydrogen atom." << endl;
-                else                 
+                else
                   *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
               }
@@ -3727,5 +3778,5 @@
           *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
         }
-  			if (Walker != OtherAtom) {	// if no white neighbours anymore, color it black 
+  			if (Walker != OtherAtom) {	// if no white neighbours anymore, color it black
           *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
   				ColorVertexList[Walker->nr] = black;
@@ -3734,5 +3785,5 @@
       }
   	} while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
-    *out << Verbose(2) << "Inner Looping is finished." << endl;    
+    *out << Verbose(2) << "Inner Looping is finished." << endl;
 
     // if we reset all AtomCount atoms, we have again technically O(N^2) ...
@@ -3750,8 +3801,8 @@
     }
   }
-  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl; 
+  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
 
   // copy together
-  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl; 
+  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
   FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
   RunningIndex = 0;
@@ -3824,18 +3875,18 @@
 
   NumCombinations = 1 << SetDimension;
-  
+
   // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
   // von Endstuecken (aus den Bonds) hinzugefï¿œï¿œgt werden und fï¿œï¿œr verbleibende ANOVAOrder
   // rekursiv GraphCrawler in der nï¿œï¿œchsten Ebene aufgerufen werden
-  
+
   *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
   *out << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->Root << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " <<  NumCombinations-1 << " combination(s)." << endl;
 
-  // initialised touched list (stores added atoms on this level) 
+  // initialised touched list (stores added atoms on this level)
   *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
   for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
     TouchedList[TouchedIndex] = -1;
   TouchedIndex = 0;
-  
+
   // create every possible combination of the endpieces
   *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
@@ -3845,5 +3896,5 @@
     for (int j=SetDimension;j--;)
       bits += (i & (1 << j)) >> j;
-      
+
     *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
     if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
@@ -3853,8 +3904,8 @@
         bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
         if (bit) {  // if bit is set, we add this bond partner
-        	OtherWalker = BondsSet[j]->rightatom;	// rightatom is always the one more distant, i.e. the one to add 
+        	OtherWalker = BondsSet[j]->rightatom;	// rightatom is always the one more distant, i.e. the one to add
           //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
-          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr); 
+          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
           if (TestKeySetInsert.second) {
             TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
@@ -3869,6 +3920,6 @@
         }
       }
-      
-      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore 
+
+      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
       if (SpaceLeft > 0) {
         *out << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
@@ -3898,5 +3949,5 @@
           }
           *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
-          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);  
+          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
           Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
         }
@@ -3907,5 +3958,5 @@
         *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
         for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
-          *out << (*runner) << " "; 
+          *out << (*runner) << " ";
         *out << endl;
         //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
@@ -3915,5 +3966,5 @@
         //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
       }
-      
+
       // --3-- remove all added items in this level from snake stack
       *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
@@ -3926,5 +3977,5 @@
       *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
       for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
-        *out << (*runner) << " "; 
+        *out << (*runner) << " ";
       *out << endl;
       TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
@@ -3933,5 +3984,5 @@
     }
   }
-  Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");  
+  Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");
   *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
 };
@@ -3942,5 +3993,5 @@
  * \return true - connected, false - disconnected
  * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
- */ 
+ */
 bool molecule::CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment)
 {
@@ -3948,8 +3999,8 @@
   bool BondStatus = false;
   int size;
-  
+
   *out << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl;
   *out << Verbose(2) << "Disconnected atom: ";
-  
+
   // count number of atoms in graph
   size = 0;
@@ -3997,10 +4048,10 @@
  * \param *out output stream for debugging
  * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
- * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on 
+ * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
  * \param RestrictedKeySet Restricted vertex set to use in context of molecule
  * \return number of inserted fragments
  * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
  */
-int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet) 
+int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
 {
   int SP, AtomKeyNr;
@@ -4023,12 +4074,12 @@
     FragmentSearch.BondsPerSPCount[i] = 0;
   FragmentSearch.BondsPerSPCount[0] = 1;
-  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root); 
+  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
   add(Binder, FragmentSearch.BondsPerSPList[1]);
-  
+
   // do a BFS search to fill the SP lists and label the found vertices
   // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
   // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
   // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
-  // (EdgeinSPLevel) of this tree ... 
+  // (EdgeinSPLevel) of this tree ...
   // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
   // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
@@ -4083,5 +4134,5 @@
     }
   }
-  
+
   // outputting all list for debugging
   *out << Verbose(0) << "Printing all found lists." << endl;
@@ -4092,7 +4143,7 @@
       Binder = Binder->next;
       *out << Verbose(2) << *Binder << endl;
-    } 
-  }
-  
+    }
+  }
+
   // creating fragments with the found edge sets  (may be done in reverse order, faster)
   SP = -1;  // the Root <-> Root edge must be subtracted!
@@ -4101,5 +4152,5 @@
     while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
       Binder = Binder->next;
-      SP ++; 
+      SP ++;
     }
   }
@@ -4108,12 +4159,12 @@
     // start with root (push on fragment stack)
     *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
-    FragmentSearch.FragmentSet->clear();  
+    FragmentSearch.FragmentSet->clear();
     *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
     // prepare the subset and call the generator
     BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
     BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond
-    
+
     SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
-    
+
     Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
   } else {
@@ -4124,5 +4175,5 @@
   // remove root from stack
   *out << Verbose(0) << "Removing root again from stack." << endl;
-  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);    
+  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
 
   // free'ing the bonds lists
@@ -4143,5 +4194,5 @@
   }
 
-  // return list  
+  // return list
   *out << Verbose(0) << "End of PowerSetGenerator." << endl;
   return (FragmentSearch.FragmentCounter - Counter);
@@ -4174,5 +4225,5 @@
     // remove bonds that are beyond bonddistance
     for(int i=NDIM;i--;)
-      Translationvector.x[i] = 0.; 
+      Translationvector.x[i] = 0.;
     // scan all bonds
     Binder = first;
@@ -4221,5 +4272,5 @@
         }
       }
-      // re-add bond    
+      // re-add bond
       link(Binder, OtherBinder);
     } else {
@@ -4275,5 +4326,5 @@
         IteratorB++;
       } // end of while loop
-    }// end of check in case of equal sizes 
+    }// end of check in case of equal sizes
   }
   return false; // if we reach this point, they are equal
@@ -4319,5 +4370,5 @@
  * \param graph1 first (dest) graph
  * \param graph2 second (source) graph
- * \param *counter keyset counter that gets increased 
+ * \param *counter keyset counter that gets increased
  */
 inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
@@ -4364,5 +4415,5 @@
   int RootKeyNr, RootNr;
   struct UniqueFragments FragmentSearch;
-  
+
   *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
 
@@ -4387,5 +4438,5 @@
     Walker = Walker->next;
     CompleteMolecule.insert(Walker->GetTrueFather()->nr);
-  } 
+  }
 
   // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
@@ -4401,10 +4452,10 @@
     //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
     //  *out << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
-    //} else 
+    //} else
     {
       // increase adaptive order by one
       Walker->GetTrueFather()->AdaptiveOrder++;
       Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
-  
+
       // initialise Order-dependent entries of UniqueFragments structure
       FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
@@ -4413,9 +4464,9 @@
         FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
         FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
-        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two 
+        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
         FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
         FragmentSearch.BondsPerSPCount[i] = 0;
-      }  
-  
+      }
+
       // allocate memory for all lower level orders in this 1D-array of ptrs
       NumLevels = 1 << (Order-1); // (int)pow(2,Order);
@@ -4423,9 +4474,9 @@
       for (int i=0;i<NumLevels;i++)
         FragmentLowerOrdersList[RootNr][i] = NULL;
-      
+
       // create top order where nothing is reduced
       *out << Verbose(0) << "==============================================================================================================" << endl;
       *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl; // , NumLevels is " << NumLevels << "
-  
+
       // Create list of Graphs of current Bond Order (i.e. F_{ij})
       FragmentLowerOrdersList[RootNr][0] =  new Graph;
@@ -4440,5 +4491,5 @@
         // we don't have to dive into suborders! These keysets are all already created on lower orders!
         // this was all ancient stuff, when we still depended on the TEFactors (and for those the suborders were needed)
-        
+
 //        if ((NumLevels >> 1) > 0) {
 //          // create lower order fragments
@@ -4447,5 +4498,5 @@
 //          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
 //            // step down to next order at (virtual) boundary of powers of 2 in array
-//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))    
+//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
 //              Order--;
 //            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
@@ -4454,5 +4505,5 @@
 //              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
 //              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
-//        
+//
 //              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
 //              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
@@ -4485,5 +4536,5 @@
       RootStack.push_back(RootKeyNr); // put back on stack
       RootNr++;
-  
+
       // free Order-dependent entries of UniqueFragments structure for next loop cycle
       Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
@@ -4491,5 +4542,5 @@
         delete(FragmentSearch.BondsPerSPList[2*i]);
         delete(FragmentSearch.BondsPerSPList[2*i+1]);
-      }  
+      }
       Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
     }
@@ -4502,6 +4553,6 @@
   Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
   delete(FragmentSearch.FragmentSet);
-  
-  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein) 
+
+  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
   // 5433222211111111
   // 43221111
@@ -4523,5 +4574,5 @@
     RootKeyNr = RootStack.front();
     RootStack.pop_front();
-    Walker = FindAtom(RootKeyNr); 
+    Walker = FindAtom(RootKeyNr);
     NumLevels = 1 << (Walker->AdaptiveOrder - 1);
     for(int i=0;i<NumLevels;i++) {
@@ -4536,5 +4587,5 @@
   Free((void **)&FragmentLowerOrdersList, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
   Free((void **)&NumMoleculesOfOrder, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
-  
+
   *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
 };
@@ -4570,6 +4621,6 @@
   atom *Walker = NULL;
   bool result = true; // status of comparison
-  
-  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl; 
+
+  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
   /// first count both their atoms and elements and update lists thereby ...
   //*out << Verbose(0) << "Counting atoms, updating list" << endl;
@@ -4618,8 +4669,8 @@
     if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
       *out << Verbose(4) << "Centers of gravity don't match." << endl;
-      result = false;  
-    }
-  }
-  
+      result = false;
+    }
+  }
+
   /// ... then make a list with the euclidian distance to this center for each atom of both molecules
   if (result) {
@@ -4637,5 +4688,5 @@
       OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
     }
-  
+
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
     *out << Verbose(5) << "Sorting distances" << endl;
@@ -4648,5 +4699,5 @@
     for(int i=AtomCount;i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
-    
+
     /// ... and compare them step by step, whether the difference is individiually(!) below \a threshold for all
     *out << Verbose(4) << "Comparing distances" << endl;
@@ -4659,5 +4710,5 @@
     Free((void **)&PermMap, "molecule::IsEqualToWithinThreshold: *PermMap");
     Free((void **)&OtherPermMap, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
-      
+
     /// free memory
     Free((void **)&Distances, "molecule::IsEqualToWithinThreshold: Distances");
@@ -4667,5 +4718,5 @@
       result = false;
     }
-  } 
+  }
   /// return pointer to map if all distances were below \a threshold
   *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
@@ -4676,5 +4727,5 @@
     *out << Verbose(3) << "Result: Not equal." << endl;
     return NULL;
-  } 
+  }
 };
 
@@ -4731,5 +4782,5 @@
  * \param *output output stream of temperature file
  * \return file written (true), failure on writing file (false)
- */ 
+ */
 bool molecule::OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output)
 {
@@ -4739,5 +4790,5 @@
   if (output == NULL)
     return false;
-  else 
+  else
     *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
   for (int step=startstep;step < endstep; step++) { // loop over all time steps
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision 1ffa216657b5beb85210757a5d292a8b64bd5b79)
+++ src/molecules.hpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -1,5 +1,5 @@
 /** \file molecules.hpp
  *
- * Class definitions of atom and molecule, element and periodentafel 
+ * Class definitions of atom and molecule, element and periodentafel
  */
 
@@ -54,15 +54,15 @@
 #define BoundariesTestPair pair< Boundaries::iterator, bool>
 
-#define PointMap map < int, class BoundaryPointSet * > 
-#define PointPair pair < int, class BoundaryPointSet * > 
-#define PointTestPair pair < PointMap::iterator, bool > 
-
-#define LineMap map < int, class BoundaryLineSet * > 
-#define LinePair pair < int, class BoundaryLineSet * > 
-#define LineTestPair pair < LinePair::iterator, bool > 
-
-#define TriangleMap map < int, class BoundaryTriangleSet * > 
-#define TrianglePair pair < int, class BoundaryTriangleSet * > 
-#define TriangleTestPair pair < TrianglePair::iterator, bool > 
+#define PointMap map < int, class BoundaryPointSet * >
+#define PointPair pair < int, class BoundaryPointSet * >
+#define PointTestPair pair < PointMap::iterator, bool >
+
+#define LineMap map < int, class BoundaryLineSet * >
+#define LinePair pair < int, class BoundaryLineSet * >
+#define LineTestPair pair < LinePair::iterator, bool >
+
+#define TriangleMap map < int, class BoundaryTriangleSet * >
+#define TrianglePair pair < int, class BoundaryTriangleSet * >
+#define TriangleTestPair pair < TrianglePair::iterator, bool >
 
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
@@ -86,5 +86,5 @@
 //bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)
 inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
-inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph 
+inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph
 int CompareDoubles (const void * a, const void * b);
 
@@ -140,13 +140,13 @@
     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();
-  
+
   bool Output(int ElementNo, int AtomNo, ofstream *out) const;
   bool OutputXYZLine(ofstream *out) const;
   atom *GetTrueFather();
   bool Compare(atom &ptr);
-  
+
   private:
 };
@@ -169,10 +169,10 @@
     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 
-  	
+    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();
@@ -180,5 +180,5 @@
   bool Contains(const atom *ptr);
   bool Contains(const int nr);
-  
+
   bond();
   bond(atom *left, atom *right);
@@ -186,6 +186,6 @@
   bond(atom *left, atom *right, int degree, int number);
   ~bond();
-    
-  private: 
+
+  private:
     enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
 };
@@ -218,8 +218,8 @@
     int NoCyclicBonds;  //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
     double BondDistance;  //!< typical bond distance used in CreateAdjacencyList() and furtheron
-  
+
   molecule(periodentafel *teil);
   ~molecule();
-  
+
   /// remove atoms from molecule.
   bool AddAtom(atom *pointer);
@@ -230,11 +230,11 @@
   atom * AddCopyAtom(atom *pointer);
   bool AddXYZFile(string filename);
-  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);  
+  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
   bond * AddBond(atom *first, atom *second, int degree);
   bool RemoveBond(bond *pointer);
   bool RemoveBonds(atom *BondPartner);
-    
+
   /// Find atoms.
-  atom * FindAtom(int Nr) const; 
+  atom * FindAtom(int Nr) const;
   atom * AskAtom(string text);
 
@@ -244,6 +244,6 @@
   void CalculateOrbitals(class config &configuration);
   bool CenterInBox(ofstream *out, Vector *BoxLengths);
-  void CenterEdge(ofstream *out, Vector *max); 
-  void CenterOrigin(ofstream *out, Vector *max); 
+  void CenterEdge(ofstream *out, Vector *max);
+  void CenterOrigin(ofstream *out, Vector *max);
   void CenterGravity(ofstream *out, Vector *max);
   void Translate(const Vector *x);
@@ -260,12 +260,13 @@
 	double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
 	bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
-	
+
   bool CheckBounds(const Vector *x) const;
   void GetAlignvector(struct lsq_params * par) const;
 
-  /// Initialising routines in fragmentation  
+  /// Initialising routines in fragmentation
+  void CreateAdjacencyList2(ofstream *out, ifstream *output);
   void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
   void CreateListOfBondsPerAtom(ofstream *out);
-  
+
   // Graph analysis
   MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
@@ -283,5 +284,5 @@
 
   molecule *CopyMolecule();
-  
+
   /// Fragment molecule by two different approaches:
   int FragmentMolecule(ofstream *out, int Order, config *configuration);
@@ -305,9 +306,9 @@
   int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
   int GuesstimateFragmentCount(ofstream *out, int order);
-  	  
-  // Recognize doubly appearing molecules in a list of them   
+
+  // Recognize doubly appearing molecules in a list of them
   int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
   int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);
-	
+
   // Output routines.
   bool Output(ofstream *out);
@@ -330,5 +331,5 @@
     int NumberOfMolecules;        //!< Number of entries in \a **FragmentList and of to be returned one.
     int NumberOfTopAtoms;         //!< Number of atoms in the molecule from which all fragments originate
-    
+
   MoleculeListClass();
   MoleculeListClass(int Num, int NumAtoms);
@@ -340,5 +341,5 @@
   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
   void Output(ofstream *out);
-  
+
   private:
 };
@@ -350,5 +351,5 @@
 class MoleculeLeafClass {
   public:
-    molecule *Leaf;                   //!< molecule of this leaf 
+    molecule *Leaf;                   //!< molecule of this leaf
     //MoleculeLeafClass *UpLeaf;        //!< Leaf one level up
     //MoleculeLeafClass *DownLeaf;      //!< First leaf one level down
@@ -386,10 +387,10 @@
     bool FastParsing;
     double Deltat;
-    
+
     private:
     char *mainname;
     char *defaultpath;
     char *pseudopotpath;
-    
+
     int DoOutVis;
     int DoOutMes;
@@ -406,5 +407,5 @@
     int UseAddGramSch;
     int Seed;
-    
+
     int MaxOuterStep;
     int OutVisStep;
@@ -414,5 +415,5 @@
     int MaxPsiStep;
     double EpsWannier;
-    
+
     int MaxMinStep;
     double RelEpsTotalEnergy;
@@ -423,7 +424,7 @@
     double InitRelEpsKineticEnergy;
     int InitMaxMinGapStopStep;
-    
+
     //double BoxLength[NDIM*NDIM];
-    
+
     double ECut;
     int MaxLevel;
@@ -434,5 +435,5 @@
     int RTActualUse;
     int AddPsis;
-    
+
     double RCut;
     int StructOpt;
@@ -441,7 +442,7 @@
     int MaxTypes;
 
-  
+
   int ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical);
-  
+
   public:
   config();
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 1ffa216657b5beb85210757a5d292a8b64bd5b79)
+++ src/vector.cpp	(revision 69eb71ca95b6ddd8113c4c58130d3a649792ac26)
@@ -1,10 +1,10 @@
 /** \file vector.cpp
- * 
+ *
  * Function implementations for the class vector.
- * 
- */
- 
+ *
+ */
+
 #include "molecules.hpp"
- 
+
 
 /************************************ Functions for class vector ************************************/
@@ -31,5 +31,5 @@
   for (int i=NDIM;i--;)
     res += (x[i]-y->x[i])*(x[i]-y->x[i]);
-  return (res);  
+  return (res);
 };
 
@@ -69,5 +69,5 @@
         if (tmp < res) res = tmp;
       }
-  return (res);  
+  return (res);
 };
 
@@ -112,5 +112,5 @@
   for (int i=NDIM;i--;)
     res += x[i]*y->x[i];
-  return (res);  
+  return (res);
 };
 
@@ -121,11 +121,11 @@
  *  \param *y array to vector with which to calculate crossproduct
  *  \return \f$ x \times y \f&
- */ 
+ */
 void Vector::VectorProduct(const Vector *y)
 {
   Vector tmp;
-  tmp[0] = x[1]*y->x[2] - x[2]*y->x[1];
-  tmp[1] = x[2]*y->x[0] - x[0]*y->x[2];
-  tmp[2] = x[0]*y->x[1] - x[1]*Y->x[0];
+  tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
+  tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
+  tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
   this->CopyVector(&tmp);
 
@@ -134,5 +134,5 @@
 
 /** projects this vector onto plane defined by \a *y.
- * \param *y array to normal vector of plane
+ * \param *y normal vector of plane
  * \return \f$\langle x, y \rangle\f$
  */
@@ -141,5 +141,6 @@
   Vector tmp;
   tmp.CopyVector(y);
-  tmp.Scale(Projection(y));
+  tmp.Normalize();
+  tmp.Scale(ScalarProduct(&tmp));
   this->SubtractVector(&tmp);
 };
@@ -162,5 +163,5 @@
   for (int i=NDIM;i--;)
     res += this->x[i]*this->x[i];
-  return (sqrt(res));  
+  return (sqrt(res));
 };
 
@@ -196,5 +197,5 @@
  */
 void Vector::Init(double x1, double x2, double x3)
-{ 
+{
   x[0] = x1;
   x[1] = x2;
@@ -202,5 +203,5 @@
 };
 
-/** Calculates the angle between this and another vector. 
+/** Calculates the angle between this and another vector.
  * \param *y array to second vector
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
@@ -208,5 +209,5 @@
 double Vector::Angle(Vector *y) const
 {
-  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
+  return acos(this->ScalarProduct(y)/Norm()/y->Norm());
 };
 
@@ -256,5 +257,5 @@
 
 /** Sums two vectors \a  and \b component-wise.
- * \param a first vector 
+ * \param a first vector
  * \param b second vector
  * \return a + b
@@ -269,5 +270,5 @@
 
 /** Factors given vector \a a times \a m.
- * \param a vector 
+ * \param a vector
  * \param m factor
  * \return a + b
@@ -327,5 +328,5 @@
 };
 
-/** Translate atom by given vector. 
+/** Translate atom by given vector.
  * \param trans[] translation vector.
  */
@@ -334,5 +335,5 @@
   for (int i=NDIM;i--;)
     x[i] += trans->x[i];
-}; 
+};
 
 /** Do a matrix multiplication.
@@ -373,5 +374,5 @@
     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-  
+
     // do the matrix multiplication
     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
@@ -397,8 +398,8 @@
 {
   for(int i=NDIM;i--;)
-    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i]; 
-};
-
-/** Mirrors atom against a given plane. 
+    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
+};
+
+/** Mirrors atom against a given plane.
  * \param n[] normal vector of mirror plane.
  */
@@ -416,5 +417,5 @@
   Output((ofstream *)&cout);
   cout << endl;
-}; 
+};
 
 /** Calculates normal vector for three given vectors (being three points in space).
@@ -448,5 +449,5 @@
   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
   Normalize();
-  
+
   return true;
 };
@@ -506,5 +507,5 @@
 /** Creates this vector as one of the possible orthonormal ones to the given one.
  * Just scan how many components of given *vector are unequal to zero and
- * try to get the skp of both to be zero accordingly.  
+ * try to get the skp of both to be zero accordingly.
  * \param *vector given vector
  * \return true - success, false - failure (null vector given)
@@ -527,5 +528,5 @@
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
-        
+
   switch(Last) {
     case 3:  // threecomponent system
@@ -540,7 +541,7 @@
     case 1: // one component system
       // set sole non-zero component to 0, and one of the other zero component pendants to 1
-      x[(Components[0]+2)%NDIM] = 0.;   
-      x[(Components[0]+1)%NDIM] = 1.;   
-      x[Components[0]] = 0.;   
+      x[(Components[0]+2)%NDIM] = 0.;
+      x[(Components[0]+1)%NDIM] = 1.;
+      x[Components[0]] = 0.;
       return true;
       break;
@@ -559,7 +560,7 @@
 {
 //  cout << Verbose(3) << "For comparison: ";
-//  cout << "A " << A->Projection(this) << "\t"; 
-//  cout << "B " << B->Projection(this) << "\t"; 
-//  cout << "C " << C->Projection(this) << "\t"; 
+//  cout << "A " << A->Projection(this) << "\t";
+//  cout << "B " << B->Projection(this) << "\t";
+//  cout << "C " << C->Projection(this) << "\t";
 //  cout << endl;
   return A->Projection(this);
@@ -571,8 +572,8 @@
  * \return true if success, false if failed due to linear dependency
  */
-bool Vector::LSQdistance(Vector **vectors, int num) 
+bool Vector::LSQdistance(Vector **vectors, int num)
 {
 	int j;
-				
+
 	for (j=0;j<num;j++) {
 		cout << Verbose(1) << j << "th atom's vector: ";
@@ -583,5 +584,5 @@
   int np = 3;
 	struct LSQ_params par;
-    
+
    const gsl_multimin_fminimizer_type *T =
      gsl_multimin_fminimizer_nmsimplex;
@@ -589,47 +590,47 @@
    gsl_vector *ss, *y;
    gsl_multimin_function minex_func;
- 
+
    size_t iter = 0, i;
    int status;
    double size;
- 
+
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
    y = gsl_vector_alloc (np);
- 
+
    /* Set all step sizes to 1 */
    gsl_vector_set_all (ss, 1.0);
- 
+
    /* Starting point */
    par.vectors = vectors;
 	 par.num = num;
-	
+
 	 for (i=NDIM;i--;)
-		gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
-	 
+		gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
+
    /* Initialize method and iterate */
    minex_func.f = &LSQ;
    minex_func.n = np;
    minex_func.params = (void *)&par;
- 
+
    s = gsl_multimin_fminimizer_alloc (T, np);
    gsl_multimin_fminimizer_set (s, &minex_func, y, 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)
          {
            printf ("converged to minimum at\n");
          }
- 
+
        printf ("%5d ", (int)iter);
        for (i = 0; i < (size_t)np; i++)
@@ -640,5 +641,5 @@
      }
    while (status == GSL_CONTINUE && iter < 100);
- 
+
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
@@ -706,7 +707,7 @@
  * \param alpha first angle
  * \param beta second angle
- * \param c norm of final vector 
+ * \param c norm of final vector
  * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
- * \bug this is not yet working properly 
+ * \bug this is not yet working properly
  */
 bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
@@ -724,5 +725,5 @@
   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
     if (fabs(x1->x[1]) > MYEPSILON) {
-      flag = 1;    
+      flag = 1;
     } else if (fabs(x1->x[2]) > MYEPSILON) {
        flag = 2;
@@ -757,11 +758,11 @@
   // now comes the case system
   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
-  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2]; 
+  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
   D3 = y->x[0]/x1->x[0]*A-B1;
   cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
   if (fabs(D1) < MYEPSILON) {
-    cout << Verbose(2) << "D1 == 0!\n"; 
+    cout << Verbose(2) << "D1 == 0!\n";
     if (fabs(D2) > MYEPSILON) {
-      cout << Verbose(3) << "D2 != 0!\n"; 
+      cout << Verbose(3) << "D2 != 0!\n";
       x[2] = -D3/D2;
       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
@@ -773,11 +774,11 @@
       cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
       if (fabs(F1) < MYEPSILON) {
-        cout << Verbose(4) << "F1 == 0!\n"; 
+        cout << Verbose(4) << "F1 == 0!\n";
         cout << Verbose(4) << "Gleichungssystem linear\n";
-        x[1] = F3/(2.*F2); 
+        x[1] = F3/(2.*F2);
       } else {
         p = F2/F1;
         q = p*p - F3/F1;
-        cout << Verbose(4) << "p " << p << "\tq " << q << endl;  
+        cout << Verbose(4) << "p " << p << "\tq " << q << endl;
         if (q < 0) {
           cout << Verbose(4) << "q < 0" << endl;
@@ -800,11 +801,11 @@
     cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     if (fabs(F1) < MYEPSILON) {
-      cout << Verbose(3) << "F1 == 0!\n"; 
+      cout << Verbose(3) << "F1 == 0!\n";
       cout << Verbose(3) << "Gleichungssystem linear\n";
-      x[2] = F3/(2.*F2);     
+      x[2] = F3/(2.*F2);
     } else {
       p = F2/F1;
       q = p*p - F3/F1;
-      cout << Verbose(3) << "p " << p << "\tq " << q << endl;  
+      cout << Verbose(3) << "p " << p << "\tq " << q << endl;
       if (q < 0) {
         cout << Verbose(3) << "q < 0" << endl;
@@ -850,5 +851,5 @@
     }
     cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
-    // apply sign matrix 
+    // apply sign matrix
     for (j=NDIM;j--;)
       x[j] *= sign[j];
@@ -856,5 +857,5 @@
     ang = x2->Angle (this);
     cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
-    if (fabs(ang - cos(beta)) < MYEPSILON) {  
+    if (fabs(ang - cos(beta)) < MYEPSILON) {
       break;
     }
