Index: src/#border.cpp#
===================================================================
--- src/#border.cpp#	(revision 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ 	(revision )
@@ -1,302 +1,0 @@
-#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 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ 	(revision )
@@ -1,10 +1,0 @@
-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 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ 	(revision )
@@ -1,463 +1,0 @@
-/** \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/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ src/boundary.cpp	(revision 12298c6c425bd81e8eb06db20b3aca630c98b764)
@@ -3,4 +3,8 @@
 
 #define DEBUG 1
+#define DoTecplotOutput 0
+#define DoRaster3DOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
 
 // ======================================== Points on Boundary =================================
@@ -564,4 +568,58 @@
 }
 ;
+
+/** Creates the objects in a raster3d file (renderable with a header.r3d)
+ * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
+{
+  atom *Walker = mol->start;
+  bond *Binder = mol->first;
+  int i;
+  Vector *center = mol->DetermineCenterOfAll(out);
+  if (rasterfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+    *rasterfile << "@header.r3d" << endl;
+    *rasterfile << "# All atoms as spheres" << endl;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      *rasterfile << "2" << endl << "  "; // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Walker->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
+
+    *rasterfile << "# All bonds as vertices" << endl;
+    while (Binder->next != mol->last) {
+      Binder = Binder->next;
+      *rasterfile << "3" << endl << "  "; // 2 is round-ended cylinder type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t";
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
+
+    *rasterfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *rasterfile << "1" << endl << "  "; // 1 is triangle type
+      for (i=0;i<3;i++) { // print each node
+        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+        *rasterfile << "\t";
+      }
+      *rasterfile << "1. 0. 0." << endl;  // red as colour
+      *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  delete(center);
+};
 
 /*
@@ -1954,5 +2012,4 @@
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out   output stream for debugging
- * @param tecplot output stream for writing found triangles in TecPlot format
  * @param mol molecule structure with all atoms and bonds
  * @param Line current baseline to search from
@@ -1961,7 +2018,7 @@
  * @param N number of found triangles
  */
-void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
+void Tesselation::Find_next_suitable_triangle(ofstream *out,
     molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    const double& RADIUS, int N)
+    const double& RADIUS, int N, const char *tempbasename)
 {
   cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
@@ -1970,5 +2027,5 @@
   Vector Chord;
   ofstream *tempstream = NULL;
-  char filename[255];
+  char NumberName[255];
   double tmp;
   //atom* Walker;
@@ -2066,24 +2123,33 @@
       cout << "Letzter Winkel bei " << TrianglesOnBoundaryCount << " Winkel ist " << Storage[2] << endl;
 
-
-  if ((TrianglesOnBoundaryCount % 1) == 0)
-    {
-    cout << Verbose(1) << "Writing temporary non convex hull to file testEnvelope-" << TriangleFilesWritten << "d.dat.\n";
-    sprintf(filename, "testEnvelope-%d.dat", TriangleFilesWritten);
-    tempstream = new ofstream(filename, ios::trunc);
-    write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten++);
-    tempstream->close();
-    tempstream->flush();
-    delete(tempstream);
+  if ((TrianglesOnBoundaryCount % 10) == 0) {
+    sprintf(NumberName, "-%d", TriangleFilesWritten);
+    if (DoTecplotOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(TecplotSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoRaster3DOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(Raster3DSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_raster3d_file(out, tempstream, this, mol);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoTecplotOutput || DoRaster3DOutput)
+      TriangleFilesWritten++;
   }
 
-  if (TrianglesOnBoundaryCount >189 and TrianglesOnBoundaryCount < 200 )
-    {
-      cerr << Verbose(0)
-          << "WARNING: Already more than " << TrianglesOnBoundaryCount-1 << "triangles, construction probably has crashed." << endl;
-      write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
-      cout << Verbose(2) << "This is currently added candidate" << Opt_Candidate << endl;
-    }
-  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
+      // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
 
   cout << Verbose(2) << " Optimal candidate is " << *Opt_Candidate << endl;
@@ -2309,5 +2375,5 @@
 ;
 
-void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
+void Find_non_convex_border(ofstream *out, const char *filename, molecule* mol)
 {
   int N = 0;
@@ -2319,5 +2385,4 @@
   cout << Verbose(0) << "Begin of Find_non_convex_border\n";
   bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
-
   Tess->Find_starting_triangle(mol, RADIUS);
 
@@ -2327,7 +2392,7 @@
       if (baseline->second->TrianglesCount == 1)
         {
-          Tess->Find_next_suitable_triangle(out, tecplot, mol,
+          Tess->Find_next_suitable_triangle(out, mol,
               *(baseline->second),
-              *(((baseline->second->triangles.begin()))->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
+              *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
           flag = true;
         }
@@ -2346,5 +2411,18 @@
     }
   cout << Verbose(1) << "Writing final tecplot file\n";
-  write_tecplot_file(out, tecplot, Tess, mol, -1);
+  if (DoTecplotOutput) {
+    string Name(filename);
+    Name.append(TecplotSuffix);
+    ofstream tecplot(Name.c_str(), ios::trunc);
+    write_tecplot_file(out, &tecplot, Tess, mol, -1);
+    tecplot.close();
+  }
+  if (DoRaster3DOutput) {
+    string Name(filename);
+    Name.append(Raster3DSuffix);
+    ofstream raster(Name.c_str(), ios::trunc);
+    write_raster3d_file(out, &raster, Tess, mol);
+    raster.close();
+  }
 
   cout << Verbose(0) << "End of Find_non_convex_border\n";
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ src/boundary.hpp	(revision 12298c6c425bd81e8eb06db20b3aca630c98b764)
@@ -86,5 +86,5 @@
     void AddTriangleToLines();
     void Find_starting_triangle(molecule* mol, const double RADIUS);
-    void Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N);
+    void Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *tempbasename);
 
     PointMap PointsOnBoundary;
@@ -111,5 +111,5 @@
 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, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule *mol, bool problem);
-void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol);
+void Find_non_convex_border(ofstream *out, const char *filename, molecule* mol);
 
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ src/builder.cpp	(revision 12298c6c425bd81e8eb06db20b3aca630c98b764)
@@ -1025,9 +1025,6 @@
             	else {
             		cout << Verbose(0) << "Evaluating npn-convex envelope.";
-            		ofstream *output = new ofstream(argv[argptr], ios::trunc);
             		cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
-            		Find_non_convex_border((ofstream *)&cout, output, mol);
-            		output->close();
-            		delete(output);
+            		Find_non_convex_border((ofstream *)&cout, argv[argptr], mol);
             		argptr+=1;
             		}
@@ -1037,5 +1034,5 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
-                cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
+                cerr << "Not enough or invalid arguments given for storing temperature: -T <temperature file>" << endl;
               } else {
                 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
