Changes in / [8ab0407:c7b1e7]

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• World.cpp (added)
• World.hpp (added)
• analysis_bonds.cpp (modified) (1 diff)
• analysis_correlation.cpp (modified) (18 diffs)
• analysis_correlation.hpp (modified) (5 diffs)
• analyzer.cpp (modified) (2 diffs)
• atom_bondedparticle.cpp (modified) (3 diffs)
• bond.cpp (modified) (2 diffs)
• bondgraph.cpp (modified) (2 diffs)
• bondgraph.hpp (modified) (1 diff)
• boundary.cpp (modified) (11 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (76 diffs)
• config.cpp (modified) (27 diffs)
• defs.hpp (modified) (1 diff)
• ellipsoid.cpp (modified) (2 diffs)
• gslvector.cpp (modified) (7 diffs)
• gslvector.hpp (modified) (3 diffs)
• helpers.hpp (modified) (1 diff)
• linkedcell.cpp (modified) (6 diffs)
• log.cpp (modified) (1 diff)
• log.hpp (modified) (1 diff)
• molecule.cpp (modified) (15 diffs)
• molecule.hpp (modified) (1 diff)
• molecule_dynamics.cpp (modified) (3 diffs)
• molecule_fragmentation.cpp (modified) (11 diffs)
• molecule_geometry.cpp (modified) (5 diffs)
• molecule_graph.cpp (modified) (4 diffs)
• moleculelist.cpp (modified) (12 diffs)
• parser.cpp (modified) (12 diffs)
• periodentafel.cpp (modified) (1 diff)
• stackclass.hpp (modified) (5 diffs)
• tesselation.cpp (modified) (75 diffs)
• tesselation.hpp (modified) (1 diff)
• tesselationhelpers.cpp (modified) (10 diffs)
• triangleintersectionlist.cpp (added)
• triangleintersectionlist.hpp (added)
• unittests/AnalysisCorrelationToSurfaceUnitTest.cpp (modified) (7 diffs)
• unittests/AnalysisCorrelationToSurfaceUnitTest.hpp (modified) (3 diffs)
• unittests/gslvectorunittest.cpp (modified) (1 diff)
• unittests/gslvectorunittest.hpp (modified) (2 diffs)
• unittests/logunittest.cpp (modified) (1 diff)
• vector.cpp (modified) (7 diffs)
• vector.hpp (modified) (2 diffs)

src/Makefile.am

@@ -8 +8 @@
 ANALYSISHEADER = analysis_bonds.hpp analysis_correlation.hpp
 
-SOURCE = ${ANALYSISSOURCE} ${ATOMSOURCE} bond.cpp bondgraph.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp errorlogger.cpp graph.cpp helpers.cpp info.cpp leastsquaremin.cpp linkedcell.cpp log.cpp logger.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
-HEADER = ${ANALYSISHEADER} ${ATOMHEADER} bond.hpp bondgraph.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp errorlogger.hpp graph.hpp helpers.hpp info.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp log.hpp logger.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
+SOURCE = ${ANALYSISSOURCE} ${ATOMSOURCE} bond.cpp bondgraph.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp errorlogger.cpp graph.cpp helpers.cpp info.cpp leastsquaremin.cpp linkedcell.cpp log.cpp logger.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp triangleintersectionlist.cpp vector.cpp verbose.cpp World.cpp
+HEADER = ${ANALYSISHEADER} ${ATOMHEADER} bond.hpp bondgraph.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp errorlogger.hpp graph.hpp helpers.hpp info.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp log.hpp logger.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp triangleintersectionlist.cpp vector.hpp verbose.hpp World.hpp
 
 BOOST_LIB = $(BOOST_LDFLAGS) $(BOOST_MPL_LIB)

src/analysis_bonds.cpp

@@ -37 +37 @@
   }
   if (((int)Mean % 2) != 0)
-    eLog() << Verbose(1) << "Something is wrong with the bond structure, the number of bonds is not even!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Something is wrong with the bond structure, the number of bonds is not even!" << endl);
   Mean /= (double)AtomCount;
 };

src/analysis_correlation.cpp

@@ -15 +15 @@
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
+#include "triangleintersectionlist.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
+#include "World.hpp"
 
 
@@ -34 +36 @@
 
   if (molecules->ListOfMolecules.empty()) {
-    eLog() << Verbose(1) <<"No molecule given." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
     return outmap;
   }
@@ -40 +42 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       atom *Walker = (*MolWalker)->start;
       while (Walker->next != (*MolWalker)->end) {
@@ -55 +57 @@
                 if (Walker->nr < OtherWalker->nr)
                   if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
+                    distance = Walker->node->PeriodicDistance(OtherWalker->node, World::get()->cell_size);
                     //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
@@ -90 +92 @@
 
   if (molecules->ListOfMolecules.empty()) {
-    eLog() << Verbose(1) <<"No molecule given." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
     return outmap;
   }
@@ -96 +98 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
-      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       atom *Walker = (*MolWalker)->start;
       while (Walker->next != (*MolWalker)->end) {
@@ -174 +176 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
+          distance = Walker->node->PeriodicDistance(point, World::get()->cell_size);
           Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
@@ -208 +210 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
@@ -255 +257 @@
 
   if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
-    Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
     return outmap;
   }
@@ -261 +263 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+      Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
-          if (triangle != NULL) {
-            distance = DistanceToTrianglePlane(Walker->node, triangle);
-            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
-          }
-        }
-      }
-    }
+          TriangleIntersectionList Intersections(Walker->node,Surface,LC);
+          distance = Intersections.GetSmallestDistance();
+          triangle = Intersections.GetClosestTriangle();
+          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
+        }
+      }
+    } else
+      Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl;
+
 
   return outmap;
@@ -312 +315 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
@@ -330 +333 @@
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
-                triangle = Surface->FindClosestTriangleToVector(&checkX, LC);
-                distance = Surface->GetDistanceSquaredToTriangle(checkX, triangle);
+                TriangleIntersectionList Intersections(&checkX,Surface,LC);
+                distance = Intersections.GetSmallestDistance();
+                triangle = Intersections.GetClosestTriangle();
                 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
                   ShortestDistance = distance;
@@ -338 +342 @@
               }
           // insert
-          ShortestDistance = sqrt(ShortestDistance);
           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
           //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
@@ -350 +353 @@
 };
 
-/** Returns the start of the bin for a given value.
+/** Returns the index of the bin for a given value.
  * \param value value whose bin to look for
  * \param BinWidth width of bin
  * \param BinStart first bin
  */
-double GetBin ( const double value, const double BinWidth, const double BinStart )
-{
-  Info FunctionInfo(__func__);
-  double bin =(double) (floor((value - BinStart)/BinWidth));
-  return (bin*BinWidth+BinStart);
+int GetBin ( const double value, const double BinWidth, const double BinStart )
+{
+  Info FunctionInfo(__func__);
+  int bin =(int) (floor((value - BinStart)/BinWidth));
+  return (bin);
 };
 
@@ -372 +375 @@
   *file << "BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << runner->first << "\t" << runner->second << endl;
+    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
   }
 };
@@ -385 +388 @@
   *file << "BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
   }
 };
@@ -400 +403 @@
     *file << runner->first;
     for (int i=0;i<NDIM;i++)
-      *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
+      *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
     *file << endl;
   }
@@ -413 +416 @@
   Info FunctionInfo(__func__);
   *file << "BinStart\tTriangle" << endl;
-  for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
-  }
-};
-
+  if (!map->empty())
+    for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
+      *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+    }
+};
+

src/analysis_correlation.hpp

@@ -51 +51 @@
 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] );
 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
-double GetBin ( const double value, const double BinWidth, const double BinStart );
+int GetBin ( const double value, const double BinWidth, const double BinStart );
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map );
 void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map );
@@ -71 +71 @@
 
   if (map == NULL) {
-    eLog() << Verbose(0) << "Nothing to min/max, map is NULL!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to min/max, map is NULL!" << endl);
     performCriticalExit();
     return;
@@ -103 +103 @@
 {
   BinPairMap *outmap = new BinPairMap;
-  double bin = 0.;
+  int bin = 0;
   double start = 0.;
   double end = 0.;
@@ -109 +109 @@
 
   if (map == NULL) {
-    eLog() << Verbose(0) << "Nothing to bin, is NULL!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to bin, is NULL!" << endl);
     performCriticalExit();
     return outmap;
@@ -122 +122 @@
     start = BinStart;
     end = BinEnd;
-    for (double runner = start; runner <= end; runner += BinWidth)
-      outmap->insert( pair<double, int> (runner, 0) );
   }
+  for (int runner = 0; runner <= ceil((end-start)/BinWidth); runner++)
+    outmap->insert( pair<double, int> ((double)runner*BinWidth+start, 0) );
 
   for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
     bin = GetBin (runner->first, BinWidth, start);
-    BinPairMapInserter = outmap->insert ( pair<double, int> (bin, 1) );
+    BinPairMapInserter = outmap->insert ( pair<double, int> ((double)bin*BinWidth+start, 1) );
     if (!BinPairMapInserter.second) {  // bin already present, increase
       BinPairMapInserter.first->second += 1;

src/analyzer.cpp

@@ -96 +96 @@
   if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
     NoHCorrection = true;
-    eLog() << Verbose(2) << "No HCorrection file found, skipping these." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "No HCorrection file found, skipping these." << endl);
   }
   
@@ -102 +102 @@
   if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
     NoHessian = true;
-    eLog() << Verbose(2) << "No Hessian file found, skipping these." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "No Hessian file found, skipping these." << endl);
   }
   if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
     NoTime = true;
-    eLog() << Verbose(2) << "No speed file found, skipping these." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "No speed file found, skipping these." << endl);
   }
   if (periode != NULL) { // also look for PAS values

src/atom_bondedparticle.cpp

@@ -86 +86 @@
       status = true;
     } else {
-      eLog() << Verbose(1) << *Binder << " does not contain " << *this << "." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << *Binder << " does not contain " << *this << "." << endl);
     }
   } else {
-    eLog() << Verbose(1) << "Binder is " << Binder << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Binder is " << Binder << "." << endl);
   }
   return status;
@@ -105 +105 @@
       status = true;
     } else {
-      eLog() << Verbose(1) << *Binder << " does not contain " << *this << "." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << *Binder << " does not contain " << *this << "." << endl);
     }
   } else {
-    eLog() << Verbose(1) << "Binder is " << Binder << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Binder is " << Binder << "." << endl);
   }
   return status;
@@ -150 +150 @@
       //Log() << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
     } else {
-      eLog() << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
+      DoeLog(2) && (eLog()<< Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl);
       FalseBondDegree++;
     }

src/bond.cpp

@@ -63 +63 @@
   if(rightatom == Atom) 
     return leftatom;
-  eLog() << Verbose(1) << "Bond " << *this << " does not contain atom " << *Atom << "!" << endl;
+  DoeLog(1) && (eLog()<< Verbose(1) << "Bond " << *this << " does not contain atom " << *Atom << "!" << endl);
   return NULL;
 };
@@ -99 +99 @@
 bool bond::MarkUsed(const enum Shading color) {
   if (Used == black) {
-    eLog() << Verbose(1) << "Bond " << this << " was already marked black!." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Bond " << this << " was already marked black!." << endl);
     return false;
   } else {

src/bondgraph.cpp

@@ -58 +58 @@
     Log() << Verbose(1) << "Parsing bond length matrix successful." << endl;
   } else {
-    eLog() << Verbose(1) << "Parsing bond length matrix failed." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Parsing bond length matrix failed." << endl);
   }
 
@@ -159 +159 @@
 {
   if (BondLengthMatrix == NULL) {// safety measure if no matrix has been parsed yet
-    eLog() << Verbose(2) << "BondLengthMatrixMinMaxDistance() called without having parsed the bond length matrix yet!" << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "BondLengthMatrixMinMaxDistance() called without having parsed the bond length matrix yet!" << endl);
     CovalentMinMaxDistance(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
   } else {

src/bondgraph.hpp

@@ -19 +19 @@
 
 #include <iostream>
+
+/*********************************************** defines ***********************************/
+
+#define BONDTHRESHOLD 0.4   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii
 
 /****************************************** forward declarations *****************************/

src/boundary.cpp

@@ -17 +17 @@
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
+#include "World.hpp"
 
 #include<gsl/gsl_poly.h>
@@ -341 +342 @@
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
         if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
-          eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl;
+          DoeLog(2) && (eLog()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl);
 
   Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
@@ -361 +362 @@
   // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
   if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
-    eLog() << Verbose(1) << "Insertion of straddling points failed!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Insertion of straddling points failed!" << endl);
 
   Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
@@ -400 +401 @@
         // flip the line
         if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
-          eLog() << Verbose(1) << "Correction of concave baselines failed!" << endl;
+          DoeLog(1) && (eLog()<< Verbose(1) << "Correction of concave baselines failed!" << endl);
         else {
           TesselStruct->FlipBaseline(line);
@@ -455 +456 @@
 
   if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
-    eLog() << Verbose(1) << "TesselStruct is empty." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty." << endl);
     return false;
   }
@@ -520 +521 @@
   // check whether there is something to work on
   if (TesselStruct == NULL) {
-    eLog() << Verbose(1) << "TesselStruct is empty!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty!" << endl);
     return volume;
   }
@@ -747 +748 @@
   Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
   if (minimumvolume > cellvolume) {
-    eLog() << Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl);
     Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
     for (int i = 0; i < NDIM; i++)
@@ -789 +790 @@
  * \param *filler molecule which the box is to be filled with
  * \param configuration contains box dimensions
+ * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
  * \param distance[NDIM] distance between filling molecules in each direction
  * \param boundary length of boundary zone between molecule and filling mollecules
@@ -797 +799 @@
  * \return *mol pointer to new molecule with filled atoms
  */
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
 {
         Info FunctionInfo(__func__);
@@ -804 +806 @@
   int N[NDIM];
   int n[NDIM];
-  double *M =  ReturnFullMatrixforSymmetric(filler->cell_size);
+  double *M =  ReturnFullMatrixforSymmetric(World::get()->cell_size);
   double Rotations[NDIM*NDIM];
+  double *MInverse = InverseMatrix(M);
   Vector AtomTranslations;
   Vector FillerTranslations;
   Vector FillerDistance;
+  Vector Inserter;
   double FillIt = false;
   atom *Walker = NULL;
   bond *Binder = NULL;
-  int i = 0;
-  LinkedCell *LCList[List->ListOfMolecules.size()];
   double phi[NDIM];
-  class Tesselation *TesselStruct[List->ListOfMolecules.size()];
-
-  i=0;
-  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
-    Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
-    LCList[i] = new LinkedCell((*ListRunner), 10.); // get linked cell list
-    Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
-    TesselStruct[i] = NULL;
-    FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
-    i++;
-  }
+  map<molecule *, Tesselation *> TesselStruct;
+  map<molecule *, LinkedCell *> LCList;
+
+  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
+    if ((*ListRunner)->AtomCount > 0) {
+      Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
+      LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
+      Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
+      TesselStruct[(*ListRunner)] = NULL;
+      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
+    }
 
   // Center filler at origin
-  filler->CenterOrigin();
+  filler->CenterEdge(&Inserter);
   filler->Center.Zero();
+  Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl;
+  Binder = filler->first;
+  while(Binder->next != filler->last) {
+    Binder = Binder->next;
+    Log() << Verbose(2) << "  " << *Binder << endl;
+  }
 
   filler->CountAtoms();
@@ -851 +859 @@
         CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
         CurrentPosition.MatrixMultiplication(M);
-        Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
-        // Check whether point is in- or outside
-        FillIt = true;
-        i=0;
-        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
-          // get linked cell list
-          if (TesselStruct[i] == NULL) {
-            eLog() << Verbose(0) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
-            FillIt = false;
+        // create molecule random translation vector ...
+        for (int i=0;i<NDIM;i++)
+          FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+        Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl;
+
+        // go through all atoms
+        for (int i=0;i<filler->AtomCount;i++)
+          CopyAtoms[i] = NULL;
+        Walker = filler->start;
+        while (Walker->next != filler->end) {
+          Walker = Walker->next;
+
+          // create atomic random translation vector ...
+          for (int i=0;i<NDIM;i++)
+            AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+
+          // ... and rotation matrix
+          if (DoRandomRotation) {
+            for (int i=0;i<NDIM;i++) {
+              phi[i] = rand()/(RAND_MAX/(2.*M_PI));
+            }
+
+            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
+            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
+            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
+            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
+            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
+            Rotations[7] =               sin(phi[1])                                                ;
+            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
+            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
+            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+          }
+
+          // ... and put at new position
+          Inserter.CopyVector(&(Walker->x));
+          if (DoRandomRotation)
+            Inserter.MatrixMultiplication(Rotations);
+          Inserter.AddVector(&AtomTranslations);
+          Inserter.AddVector(&FillerTranslations);
+          Inserter.AddVector(&CurrentPosition);
+
+          // check whether inserter is inside box
+          Inserter.MatrixMultiplication(MInverse);
+          FillIt = true;
+          for (int i=0;i<NDIM;i++)
+            FillIt = FillIt && (Inserter.x[i] >= -MYEPSILON) && ((Inserter.x[i]-1.) <= MYEPSILON);
+          Inserter.MatrixMultiplication(M);
+
+          // Check whether point is in- or outside
+          for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
+            // get linked cell list
+            if (TesselStruct[(*ListRunner)] != NULL) {
+              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
+              FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
+            }
+          }
+          // insert into Filling
+          if (FillIt) {
+            Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl;
+            // copy atom ...
+            CopyAtoms[Walker->nr] = new atom(Walker);
+            CopyAtoms[Walker->nr]->x.CopyVector(&Inserter);
+            Filling->AddAtom(CopyAtoms[Walker->nr]);
+            Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
           } else {
-            const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
-            FillIt = FillIt && (distance > boundary*boundary);
-            if (FillIt) {
-              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
-            } else {
-              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
-              break;
-            }
-            i++;
+            Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl;
+            CopyAtoms[Walker->nr] = NULL;
+            continue;
           }
         }
-
-        if (FillIt) {
-          // fill in Filler
-          Log() << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
-
-          // create molecule random translation vector ...
-          for (int i=0;i<NDIM;i++)
-            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-          Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
-
-          // go through all atoms
-          Walker = filler->start;
-          while (Walker->next != filler->end) {
-            Walker = Walker->next;
-            // copy atom ...
-            CopyAtoms[Walker->nr] = new atom(Walker);
-
-            // create atomic random translation vector ...
-            for (int i=0;i<NDIM;i++)
-              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-
-            // ... and rotation matrix
-            if (DoRandomRotation) {
-              for (int i=0;i<NDIM;i++) {
-                phi[i] = rand()/(RAND_MAX/(2.*M_PI));
-              }
-
-              Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
-              Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
-              Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
-              Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
-              Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
-              Rotations[7] =               sin(phi[1])                                                ;
-              Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
-              Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
-              Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
-            }
-
-            // ... and put at new position
-            if (DoRandomRotation)
-              CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
-            CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
-            CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
-            CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
-
-            // insert into Filling
-
-            // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
-            Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
-            Filling->AddAtom(CopyAtoms[Walker->nr]);
-          }
-
-          // go through all bonds and add as well
-          Binder = filler->first;
-          while(Binder->next != filler->last) {
-            Binder = Binder->next;
-            Log() << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
+        // go through all bonds and add as well
+        Binder = filler->first;
+        while(Binder->next != filler->last) {
+          Binder = Binder->next;
+          if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
+            Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
             Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
           }
-        } else {
-          // leave empty
-          Log() << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
         }
       }
   Free(&M);
-
-  // output to file
-  TesselStruct[0]->LastTriangle = NULL;
-  StoreTrianglesinFile(Filling, TesselStruct[0], "Tesselated", ".dat");
-
-  for (size_t i=0;i<List->ListOfMolecules.size();i++) {
-        delete(LCList[i]);
-        delete(TesselStruct[i]);
-  }
+  Free(&MInverse);
+
   return Filling;
 };

src/boundary.hpp

@@ -49 +49 @@
 
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
 void FindConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);

src/builder.cpp

@@ -68 +68 @@
 #include "periodentafel.hpp"
 #include "version.h"
+#include "World.hpp"
 
 /********************************************* Subsubmenu routine ************************************/
@@ -84 +85 @@
   bool valid;
 
-  Log() << Verbose(0) << "===========ADD ATOM============================" << endl;
-  Log() << Verbose(0) << " a - state absolute coordinates of atom" << endl;
-  Log() << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
-  Log() << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
-  Log() << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
-  Log() << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
-  Log() << Verbose(0) << "all else - go back" << endl;
-  Log() << Verbose(0) << "===============================================" << endl;
-  Log() << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
-  Log() << Verbose(0) << "INPUT: ";
+  cout << Verbose(0) << "===========ADD ATOM============================" << endl;
+  cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
+  cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
+  cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
+  cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
+  cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
+  cout << Verbose(0) << "INPUT: ";
   cin >> choice;
 
   switch (choice) {
     default:
-      eLog() << Verbose(2) << "Not a valid choice." << endl;
+      DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl);
       break;
       case 'a': // absolute coordinates of atom
-        Log() << Verbose(0) << "Enter absolute coordinates." << endl;
+        cout << Verbose(0) << "Enter absolute coordinates." << endl;
         first = new atom;
-        first->x.AskPosition(mol->cell_size, false);
+        first->x.AskPosition(World::get()->cell_size, false);
         first->type = periode->AskElement();  // give type
         mol->AddAtom(first);  // add to molecule
@@ -112 +113 @@
         valid = true;
         do {
-          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
-          Log() << Verbose(0) << "Enter reference coordinates." << endl;
-          x.AskPosition(mol->cell_size, true);
-          Log() << Verbose(0) << "Enter relative coordinates." << endl;
-          first->x.AskPosition(mol->cell_size, false);
+          if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
+          cout << Verbose(0) << "Enter reference coordinates." << endl;
+          x.AskPosition(World::get()->cell_size, true);
+          cout << Verbose(0) << "Enter relative coordinates." << endl;
+          first->x.AskPosition(World::get()->cell_size, false);
           first->x.AddVector((const Vector *)&x);
-          Log() << Verbose(0) << "\n";
+          cout << Verbose(0) << "\n";
         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
         first->type = periode->AskElement();  // give type
@@ -128 +129 @@
         valid = true;
         do {
-          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
+          if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
           second = mol->AskAtom("Enter atom number: ");
           Log() << Verbose(0) << "Enter relative coordinates." << endl;
-          first->x.AskPosition(mol->cell_size, false);
+          first->x.AskPosition(World::get()->cell_size, false);
           for (int i=NDIM;i--;) {
             first->x.x[i] += second->x.x[i];
@@ -145 +146 @@
         do {
           if (!valid) {
-            eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl;
+            DoeLog(2) && (eLog()<< Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl);
           }
-          Log() << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
+          cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
           second = mol->AskAtom("Enter central atom: ");
           third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
@@ -158 +159 @@
           c *= M_PI/180.;
           bound(&c, -M_PI, M_PI);
-          Log() << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
+          cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
 /*
           second->Output(1,1,(ofstream *)&cout);
@@ -170 +171 @@
 
           if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
-            Log() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
+         coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
             continue;
           }
@@ -247 +248 @@
           atoms[i] = NULL;
         int i=0, j=0;
-        Log() << Verbose(0) << "Now we need at least three molecules.\n";
+        cout << Verbose(0) << "Now we need at least three molecules.\n";
         do {
-          Log() << Verbose(0) << "Enter " << i+1 << "th atom: ";
+          cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
           cin >> j;
           if (j != -1) {
@@ -264 +265 @@
         } else {
           delete first;
-          Log() << Verbose(0) << "Please enter at least two vectors!\n";
+          cout << Verbose(0) << "Please enter at least two vectors!\n";
         }
         break;
@@ -278 +279 @@
   char choice;  // menu choice char
 
-  Log() << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
-  Log() << Verbose(0) << " a - on origin" << endl;
-  Log() << Verbose(0) << " b - on center of gravity" << endl;
-  Log() << Verbose(0) << " c - within box with additional boundary" << endl;
-  Log() << Verbose(0) << " d - within given simulation box" << endl;
-  Log() << Verbose(0) << "all else - go back" << endl;
-  Log() << Verbose(0) << "===============================================" << endl;
-  Log() << Verbose(0) << "INPUT: ";
+  cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - on origin" << endl;
+  cout << Verbose(0) << " b - on center of gravity" << endl;
+  cout << Verbose(0) << " c - within box with additional boundary" << endl;
+  cout << Verbose(0) << " d - within given simulation box" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
   cin >> choice;
 
   switch (choice) {
     default:
-      Log() << Verbose(0) << "Not a valid choice." << endl;
+      cout << Verbose(0) << "Not a valid choice." << endl;
       break;
     case 'a':
-      Log() << Verbose(0) << "Centering atoms in config file on origin." << endl;
+      cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
       mol->CenterOrigin();
       break;
     case 'b':
-      Log() << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
+      cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
       mol->CenterPeriodic();
       break;
     case 'c':
-      Log() << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
+      cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
       for (int i=0;i<NDIM;i++) {
-        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
+        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
         cin >> y.x[i];
       }
@@ -314 +315 @@
       break;
     case 'd':
-      Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+      cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
       for (int i=0;i<NDIM;i++) {
-        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
+        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
         cin >> x.x[i];
       }
@@ -337 +338 @@
   char choice;  // menu choice char
 
-  Log() << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
-  Log() << Verbose(0) << " a - state three atoms defining align plane" << endl;
-  Log() << Verbose(0) << " b - state alignment vector" << endl;
-  Log() << Verbose(0) << " c - state two atoms in alignment direction" << endl;
-  Log() << Verbose(0) << " d - align automatically by least square fit" << endl;
-  Log() << Verbose(0) << "all else - go back" << endl;
-  Log() << Verbose(0) << "===============================================" << endl;
-  Log() << Verbose(0) << "INPUT: ";
+  cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
+  cout << Verbose(0) << " b - state alignment vector" << endl;
+  cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
+  cout << Verbose(0) << " d - align automatically by least square fit" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
   cin >> choice;
 
@@ -357 +358 @@
       break;
     case 'b': // normal vector of mirror plane
-      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
-      n.AskPosition(mol->cell_size,0);
+      cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
+      n.AskPosition(World::get()->cell_size,0);
       n.Normalize();
       break;
@@ -377 +378 @@
         fscanf(stdin, "%3s", shorthand);
       } while ((param.type = periode->FindElement(shorthand)) == NULL);
-      Log() << Verbose(0) << "Element is " << param.type->name << endl;
+      cout << Verbose(0) << "Element is " << param.type->name << endl;
       mol->GetAlignvector(&param);
       for (int i=NDIM;i--;) {
@@ -384 +385 @@
       }
       gsl_vector_free(param.x);
-      Log() << Verbose(0) << "Offset vector: ";
+      cout << Verbose(0) << "Offset vector: ";
       x.Output();
       Log() << Verbose(0) << endl;
@@ -425 +426 @@
     case 'b': // normal vector of mirror plane
       Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
-      n.AskPosition(mol->cell_size,0);
+      n.AskPosition(World::get()->cell_size,0);
       n.Normalize();
       break;
@@ -655 +656 @@
 static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
 {
-  atom *first, *second;
+  atom *first, *second, *third;
   molecule *mol = NULL;
   Vector x,y,z,n; // coordinates for absolute point in cell volume
@@ -667 +668 @@
   Log() << Verbose(0) << "r - remove an atom" << endl;
   Log() << Verbose(0) << "b - scale a bond between atoms" << endl;
+  Log() << Verbose(0) << "t - turn an atom round another bond" << endl;
   Log() << Verbose(0) << "u - change an atoms element" << endl;
   Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
@@ -672 +674 @@
   Log() << Verbose(0) << "===============================================" << endl;
   if (molecules->NumberOfActiveMolecules() > 1)
-    eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
   Log() << Verbose(0) << "INPUT: ";
   cin >> choice;
@@ -748 +750 @@
       break;
 
+    case 't': // turn/rotate atom
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
+          mol = *ListRunner;
+          Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl;
+          first = mol->AskAtom("Enter turning atom: ");
+          second = mol->AskAtom("Enter central atom: ");
+          third  = mol->AskAtom("Enter bond atom: ");
+          cout << Verbose(0) << "Enter new angle in degrees: ";
+          double tmp = 0.;
+          cin >> tmp;
+          // calculate old angle
+          x.CopyVector((const Vector *)&first->x);
+          x.SubtractVector((const Vector *)&second->x);
+          y.CopyVector((const Vector *)&third->x);
+          y.SubtractVector((const Vector *)&second->x);
+          double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
+          cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
+          cout << Verbose(0) << alpha << " degrees" << endl;
+          // rotate
+          z.MakeNormalVector(&x,&y);
+          x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
+          x.AddVector(&second->x);
+          first->x.CopyVector(&x);
+          // check new angle
+          x.CopyVector((const Vector *)&first->x);
+          x.SubtractVector((const Vector *)&second->x);
+          alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
+          cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
+          cout << Verbose(0) << alpha << " degrees" << endl;
+        }
+      break;
+
     case 'u': // change an atom's element
       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
@@ -795 +830 @@
   Log() << Verbose(0) << "===============================================" << endl;
   if (molecules->NumberOfActiveMolecules() > 1)
-    eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
   Log() << Verbose(0) << "INPUT: ";
   cin >> choice;
@@ -828 +863 @@
           }
           if (count != j)
-            eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+            DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
           x.Zero();
           y.Zero();
-          y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+          y.x[abs(axis)-1] = World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
           for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
             x.AddVector(&y); // per factor one cell width further
@@ -854 +889 @@
             mol->Translate(&x);
           }
-          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
+          World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
         }
       }
@@ -911 +946 @@
         Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
         Log() << Verbose(0) << "Enter translation vector." << endl;
-        x.AskPosition(mol->cell_size,0);
+        x.AskPosition(World::get()->cell_size,0);
         mol->Center.AddVector((const Vector *)&x);
      }
@@ -971 +1006 @@
         // center at set box dimensions
         mol->CenterEdge(&center);
-        mol->cell_size[0] = center.x[0];
-        mol->cell_size[1] = 0;
-        mol->cell_size[2] = center.x[1];
-        mol->cell_size[3] = 0;
-        mol->cell_size[4] = 0;
-        mol->cell_size[5] = center.x[2];
+        double * const cell_size = World::get()->cell_size;
+        cell_size[0] = center.x[0];
+        cell_size[1] = 0;
+        cell_size[2] = center.x[1];
+        cell_size[3] = 0;
+        cell_size[4] = 0;
+        cell_size[5] = center.x[2];
         molecules->insert(mol);
       }
@@ -1172 +1208 @@
     mol = (molecules->ListOfMolecules.front())->CopyMolecule();
   else {
-    eLog() << Verbose(0) << "I don't have anything to test on ... ";
+    DoeLog(0) && (eLog()<< Verbose(0) << "I don't have anything to test on ... ");
     performCriticalExit();
     return;
@@ -1253 +1289 @@
 
   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
-    eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
   }
 
@@ -1357 +1393 @@
 
   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
-    eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
   }
 
@@ -1387 +1423 @@
   enum ConfigStatus configPresent = absent;
   clock_t start,end;
+  double MaxDistance = -1;
   int argptr;
   molecule *mol = NULL;
@@ -1412 +1449 @@
             Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
             Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
-            Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
+            Log() << Verbose(0) << "\t-C <type> [params] <output> <bin output> <BinWidth> <BinStart> <BinEnd>\tPair Correlation analysis." << endl;
             Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
@@ -1418 +1455 @@
             Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
             Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
-            Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
+            Log() << Verbose(0) << "\t-F <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
+            Log() << Verbose(0) << "\t-FF <MaxDistance> <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
             Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
             Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
@@ -1442 +1480 @@
             Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl;
             Log() << Verbose(0) << "\t-V\t\tGives version information." << endl;
+            Log() << Verbose(0) << "\t-X\t\tset default name of a molecule." << endl;
             Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl;
             return (1);
@@ -1457 +1496 @@
             return (1);
             break;
+          case 'B':
+            if (ExitFlag == 0) ExitFlag = 1;
+            if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
+              ExitFlag = 255;
+              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
+              performCriticalExit();
+            } else {
+              SaveFlag = true;
+              j = -1;
+              Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+              double * const cell_size = World::get()->cell_size;
+              for (int i=0;i<6;i++) {
+                cell_size[i] = atof(argv[argptr+i]);
+              }
+              argptr+=6;
+            }
+            break;
           case 'e':
             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-              eLog() << Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
+              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl);
               performCriticalExit();
             } else {
@@ -1469 +1525 @@
           case 'g':
             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-              eLog() << Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl;
+              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl);
               performCriticalExit();
             } else {
@@ -1480 +1536 @@
             Log() << Verbose(0) << "I won't parse trajectories." << endl;
             configuration.FastParsing = true;
+            break;
+          case 'X':
+            {
+              char **name = &(World::get()->DefaultName);
+              delete[](*name);
+              const int length = strlen(argv[argptr]);
+              *name = new char[length+2];
+              strncpy(*name, argv[argptr], length);
+              Log() << Verbose(0) << "Default name of new molecules set to " << *name << "." << endl;
+            }
             break;
           default:   // no match? Step on
@@ -1556 +1622 @@
          Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
        } else {
-         eLog() << Verbose(1) << "Bond length table loading failed." << endl;
+         DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
        }
      }
@@ -1572 +1638 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl);
                 performCriticalExit();
               } else {
@@ -1589 +1655 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3]))) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl);
                 performCriticalExit();
               } else {
@@ -1605 +1671 @@
                     configPresent = present;
                 } else
-                  eLog() << Verbose(1) << "Could not find the specified element." << endl;
+                  DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the specified element." << endl);
                 argptr+=4;
               }
@@ -1618 +1684 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl);
                 performCriticalExit();
               } else {
@@ -1671 +1737 @@
                   }
               }
-              if (mol == NULL) {
+              if ((mol == NULL) && (!molecules->ListOfMolecules.empty())) {
                 mol = *(molecules->ListOfMolecules.begin());
-                mol->ActiveFlag = true;
+                if (mol != NULL)
+                  mol->ActiveFlag = true;
               }
               break;
             case 'C':
-              if (ExitFlag == 0) ExitFlag = 1;
-              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) {
-                ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl;
-                performCriticalExit();
-              } else {
-                ofstream output(argv[argptr+1]);
-                ofstream binoutput(argv[argptr+2]);
-                const double radius = 5.;
-
-                // get the boundary
-                class molecule *Boundary = NULL;
-                class Tesselation *TesselStruct = NULL;
-                const LinkedCell *LCList = NULL;
-                // find biggest molecule
-                int counter  = 0;
-                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
-                  if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
-                    Boundary = *BigFinder;
+              {
+                int ranges[3] = {1, 1, 1};
+                bool periodic = (argv[argptr-1][2] =='p');
+                if (ExitFlag == 0) ExitFlag = 1;
+                if ((argptr >= argc)) {
+                  ExitFlag = 255;
+                  DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C[p] <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl);
+                  performCriticalExit();
+                } else {
+                  switch(argv[argptr][0]) {
+                    case 'E':
+                      {
+                        if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) {
+                          ExitFlag = 255;
+                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl);
+                          performCriticalExit();
+                        } else {
+                          ofstream output(argv[argptr+3]);
+                          ofstream binoutput(argv[argptr+4]);
+                          const double BinStart = atof(argv[argptr+5]);
+                          const double BinEnd = atof(argv[argptr+6]);
+
+                          element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
+                          element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2]));
+                          PairCorrelationMap *correlationmap = NULL;
+                          if (periodic)
+                            correlationmap = PeriodicPairCorrelation(molecules, elemental, elemental2, ranges);
+                          else
+                            correlationmap = PairCorrelation(molecules, elemental, elemental2);
+                          //OutputCorrelationToSurface(&output, correlationmap);
+                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
+                          OutputCorrelation ( &binoutput, binmap );
+                          output.close();
+                          binoutput.close();
+                          delete(binmap);
+                          delete(correlationmap);
+                          argptr+=7;
+                        }
+                      }
+                      break;
+
+                    case 'P':
+                      {
+                        if ((argptr+8 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+7])) || (!IsValidNumber(argv[argptr+8])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) {
+                          ExitFlag = 255;
+                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl);
+                          performCriticalExit();
+                        } else {
+                          ofstream output(argv[argptr+5]);
+                          ofstream binoutput(argv[argptr+6]);
+                          const double BinStart = atof(argv[argptr+7]);
+                          const double BinEnd = atof(argv[argptr+8]);
+
+                          element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
+                          Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3]));
+                          CorrelationToPointMap *correlationmap = NULL;
+                          if (periodic)
+                            correlationmap  = PeriodicCorrelationToPoint(molecules, elemental, Point, ranges);
+                          else
+                            correlationmap = CorrelationToPoint(molecules, elemental, Point);
+                          //OutputCorrelationToSurface(&output, correlationmap);
+                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
+                          OutputCorrelation ( &binoutput, binmap );
+                          output.close();
+                          binoutput.close();
+                          delete(Point);
+                          delete(binmap);
+                          delete(correlationmap);
+                          argptr+=9;
+                        }
+                      }
+                      break;
+
+                    case 'S':
+                      {
+                        if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) {
+                          ExitFlag = 255;
+                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C S <Z> <output> <bin output> <BinWidth> <BinStart> <BinEnd>" << endl);
+                          performCriticalExit();
+                        } else {
+                          ofstream output(argv[argptr+2]);
+                          ofstream binoutput(argv[argptr+3]);
+                          const double radius = 4.;
+                          const double BinWidth = atof(argv[argptr+4]);
+                          const double BinStart = atof(argv[argptr+5]);
+                          const double BinEnd = atof(argv[argptr+6]);
+                          double LCWidth = 20.;
+                          if (BinEnd > 0) {
+                            if (BinEnd > 2.*radius)
+                                LCWidth = BinEnd;
+                            else
+                              LCWidth = 2.*radius;
+                          }
+
+                          // get the boundary
+                          class molecule *Boundary = NULL;
+                          class Tesselation *TesselStruct = NULL;
+                          const LinkedCell *LCList = NULL;
+                          // find biggest molecule
+                          int counter  = 0;
+                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
+                            if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
+                              Boundary = *BigFinder;
+                            }
+                            counter++;
+                          }
+                          bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
+                          counter = 0;
+                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
+                            Actives[counter++] = (*BigFinder)->ActiveFlag;
+                            (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
+                          }
+                          LCList = new LinkedCell(Boundary, LCWidth);
+                          element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
+                          FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
+                          CorrelationToSurfaceMap *surfacemap = NULL;
+                          if (periodic)
+                            surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges);
+                          else
+                            surfacemap = CorrelationToSurface( molecules, elemental, TesselStruct, LCList);
+                          OutputCorrelationToSurface(&output, surfacemap);
+                          // check whether radius was appropriate
+                          {
+                            double start; double end;
+                            GetMinMax( surfacemap, start, end);
+                            if (LCWidth < end)
+                              DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl);
+                          }
+                          BinPairMap *binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd );
+                          OutputCorrelation ( &binoutput, binmap );
+                          output.close();
+                          binoutput.close();
+                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
+                            (*BigFinder)->ActiveFlag = Actives[counter++];
+                          Free(&Actives);
+                          delete(LCList);
+                          delete(TesselStruct);
+                          delete(binmap);
+                          delete(surfacemap);
+                          argptr+=7;
+                        }
+                      }
+                      break;
+
+                    default:
+                      ExitFlag = 255;
+                      DoeLog(0) && (eLog()<< Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl);
+                      performCriticalExit();
+                      break;
                   }
-                  counter++;
                 }
-                bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
-                counter = 0;
-                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
-                  Actives[counter++] = (*BigFinder)->ActiveFlag;
-                  (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
-                }
-                LCList = new LinkedCell(Boundary, 2.*radius);
-                element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
-                FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
-                int ranges[NDIM] = {1,1,1};
-                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
-                OutputCorrelationToSurface(&output, surfacemap);
-                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 20. );
-                OutputCorrelation ( &binoutput, binmap );
-                output.close();
-                binoutput.close();
-                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
-                  (*BigFinder)->ActiveFlag = Actives[counter++];
-                Free(&Actives);
-                delete(LCList);
-                delete(TesselStruct);
-                argptr+=3;
-              }
-              break;
+                break;
+              }
             case 'E':
               if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl);
                 performCriticalExit();
               } else {
@@ -1739 +1913 @@
             case 'F':
               if (ExitFlag == 0) ExitFlag = 1;
-              if (argptr+6 >=argc) {
+              MaxDistance = -1;
+              if (argv[argptr-1][2] == 'F') { // option is -FF?
+                // fetch first argument as max distance to surface
+                MaxDistance = atof(argv[argptr++]);
+                Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl;
+              }
+              if ((argptr+7 >=argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+7]))) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <xyz of filler> <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl);
                 performCriticalExit();
               } else {
@@ -1748 +1928 @@
                 // construct water molecule
                 molecule *filler = new molecule(periode);
+                if (!filler->AddXYZFile(argv[argptr])) {
+                  DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse filler molecule from " << argv[argptr] << "." << endl);
+                }
+                filler->SetNameFromFilename(argv[argptr]);
+                configuration.BG->ConstructBondGraph(filler);
                 molecule *Filling = NULL;
-                atom *second = NULL, *third = NULL;
-//                first = new atom();
-//                first->type = periode->FindElement(5);
-//                first->x.Zero();
-//                filler->AddAtom(first);
-                first = new atom();
-                first->type = periode->FindElement(1);
-                first->x.Init(0.441, -0.143, 0.);
-                filler->AddAtom(first);
-                second = new atom();
-                second->type = periode->FindElement(1);
-                second->x.Init(-0.464, 1.137, 0.0);
-                filler->AddAtom(second);
-                third = new atom();
-                third->type = periode->FindElement(8);
-                third->x.Init(-0.464, 0.177, 0.);
-                filler->AddAtom(third);
-                filler->AddBond(first, third, 1);
-                filler->AddBond(second, third, 1);
                 // call routine
                 double distance[NDIM];
                 for (int i=0;i<NDIM;i++)
-                  distance[i] = atof(argv[argptr+i]);
-                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
+                  distance[i] = atof(argv[argptr+i+1]);
+                Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+4]), atof(argv[argptr+5]), atof(argv[argptr+6]), atoi(argv[argptr+7]));
                 if (Filling != NULL) {
                   Filling->ActiveFlag = false;
@@ -1785 +1951 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag =255;
-                eLog() << Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl);
                 performCriticalExit();
               } else {
@@ -1800 +1966 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag =255;
-                eLog() << Verbose(0) << "Missing path of adjacency file: -j <path>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of adjacency file: -j <path>" << endl);
                 performCriticalExit();
               } else {
                 Log() << Verbose(0) << "Storing adjacency to path " << argv[argptr] << "." << endl;
                 configuration.BG->ConstructBondGraph(mol);
-                mol->StoreAdjacencyToFile(argv[argptr]);
+                mol->StoreAdjacencyToFile(NULL, argv[argptr]);
                 argptr+=1;
               }
@@ -1814 +1980 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag =255;
-                eLog() << Verbose(0) << "Missing path of bonds file: -j <path>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of bonds file: -j <path>" << endl);
                 performCriticalExit();
               } else {
                 Log() << Verbose(0) << "Storing bonds to path " << argv[argptr] << "." << endl;
                 configuration.BG->ConstructBondGraph(mol);
-                mol->StoreBondsToFile(argv[argptr]);
+                mol->StoreBondsToFile(NULL, argv[argptr]);
                 argptr+=1;
               }
@@ -1828 +1994 @@
               if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl);
                 performCriticalExit();
               } else {
@@ -1864 +2030 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl);
                 performCriticalExit();
               } else {
@@ -1882 +2048 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl);
                 performCriticalExit();
               } else {
@@ -1900 +2066 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl);
                 performCriticalExit();
               } else {
@@ -1916 +2082 @@
               if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl);
                 performCriticalExit();
               } else {
@@ -1933 +2099 @@
                   }
                 } else {
-                  eLog() << Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl;
+                  DoeLog(1) && (eLog()<< Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl);
                 }
                 argptr+=2;
@@ -1942 +2108 @@
               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl);
                 performCriticalExit();
               } else {
@@ -1958 +2124 @@
               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl);
                 performCriticalExit();
               } else {
@@ -1974 +2140 @@
               if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl);
                 performCriticalExit();
               } else {
@@ -1985 +2151 @@
                 factor[2] = atof(argv[argptr+2]);
                 mol->Scale((const double ** const)&factor);
+                double * const cell_size = World::get()->cell_size;
                 for (int i=0;i<NDIM;i++) {
                   j += i+1;
                   x.x[i] = atof(argv[NDIM+i]);
-                  mol->cell_size[j]*=factor[i];
+                  cell_size[j]*=factor[i];
                 }
                 delete[](factor);
@@ -1998 +2165 @@
               if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
                 performCriticalExit();
               } else {
@@ -2004 +2171 @@
                 j = -1;
                 Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+                double * const cell_size = World::get()->cell_size;
                 for (int i=0;i<6;i++) {
-                  mol->cell_size[i] = atof(argv[argptr+i]);
+                  cell_size[i] = atof(argv[argptr+i]);
                 }
                 // center
@@ -2016 +2184 @@
               if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
                 performCriticalExit();
               } else {
@@ -2022 +2190 @@
                 j = -1;
                 Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+                double * const cell_size = World::get()->cell_size;
                 for (int i=0;i<6;i++) {
-                  mol->cell_size[i] = atof(argv[argptr+i]);
+                  cell_size[i] = atof(argv[argptr+i]);
                 }
                 // center
@@ -2034 +2203 @@
               if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl);
                 performCriticalExit();
               } else {
@@ -2045 +2214 @@
                 mol->SetBoxDimension(&x);
                 // translate each coordinate by boundary
+                double * const cell_size = World::get()->cell_size;
                 j=-1;
                 for (int i=0;i<NDIM;i++) {
                   j += i+1;
                   x.x[i] = atof(argv[argptr+i]);
-                  mol->cell_size[j] += x.x[i]*2.;
+                  cell_size[j] += x.x[i]*2.;
                 }
                 mol->Translate((const Vector *)&x);
@@ -2068 +2238 @@
               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])))  {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl);
                 performCriticalExit();
               } else {
@@ -2082 +2252 @@
               if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl);
                 performCriticalExit();
               } else {
@@ -2102 +2272 @@
               j = atoi(argv[argptr++]);
               if ((j<0) || (j>1)) {
-                eLog() << Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
+                DoeLog(1) && (eLog()<< Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl);
                 j = 0;
               }
@@ -2116 +2286 @@
               if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl);
                 performCriticalExit();
               } else {
@@ -2141 +2311 @@
               if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl);
                 performCriticalExit();
               } else {
@@ -2152 +2322 @@
                 if (volume != -1)
                   ExitFlag = 255;
-                  eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;
+                  DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl);
                   performCriticalExit();
               } else {
@@ -2160 +2330 @@
                 density = atof(argv[argptr++]);
                 if (density < 1.0) {
-                  eLog() << Verbose(1) << "Density must be greater than 1.0g/cm^3 !" << endl;
+                  DoeLog(1) && (eLog()<< Verbose(1) << "Density must be greater than 1.0g/cm^3 !" << endl);
                   density = 1.3;
                 }
@@ -2166 +2336 @@
 //                  repetition[i] = atoi(argv[argptr++]);
 //                  if (repetition[i] < 1)
-//                    eLog() << Verbose(1) << "repetition value must be greater 1!" << endl;
+//                    DoeLog(1) && (eLog()<< Verbose(1) << "repetition value must be greater 1!" << endl);
 //                  repetition[i] = 1;
 //                }
@@ -2176 +2346 @@
               if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl);
                 performCriticalExit();
               } else {
                 SaveFlag = true;
+                double * const cell_size = World::get()->cell_size;
                 for (int axis = 1; axis <= NDIM; axis++) {
                   int faktor = atoi(argv[argptr++]);
@@ -2186 +2357 @@
                   Vector ** vectors;
                   if (faktor < 1) {
-                    eLog() << Verbose(1) << "Repetition factor mus be greater than 1!" << endl;
+                    DoeLog(1) && (eLog()<< Verbose(1) << "Repetition factor mus be greater than 1!" << endl);
                     faktor = 1;
                   }
@@ -2203 +2374 @@
                     }
                     if (count != j)
-                      eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+                      DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
                     x.Zero();
                     y.Zero();
-                    y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+                    y.x[abs(axis)-1] = cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
                     for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
                       x.AddVector(&y); // per factor one cell width further
@@ -2227 +2398 @@
                       mol->Translate(&x);
                     }
-                    mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
+                    cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
                   }
                 }
@@ -2269 +2440 @@
 
   cout << ESPACKVersion << endl;
+
+  Log() << Verbose(1) << "test" << endl;
+  DoLog(3) && (Log() << Verbose(1) << "test");
 
   setVerbosity(0);
@@ -2296 +2470 @@
   if (molecules->ListOfMolecules.size() == 0) {
     mol = new molecule(periode);
-    if (mol->cell_size[0] == 0.) {
+    double * const cell_size = World::get()->cell_size;
+    if (cell_size[0] == 0.) {
       Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
       for (int i=0;i<6;i++) {
         Log() << Verbose(1) << "Cell size" << i << ": ";
-        cin >> mol->cell_size[i];
+        cin >> cell_size[i];
       }
     }

src/config.cpp

@@ -19 +19 @@
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "World.hpp"
 
 /******************************** Functions for class ConfigFileBuffer **********************/
@@ -73 +74 @@
   file= new ifstream(filename);
   if (file == NULL) {
-    eLog() << Verbose(1) << "config file " << filename << " missing!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
     return;
   }
@@ -88 +89 @@
   // allocate buffer's 1st dimension
   if (buffer != NULL) {
-    eLog() << Verbose(1) << "FileBuffer->buffer is not NULL!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "FileBuffer->buffer is not NULL!" << endl);
     return;
   } else
@@ -143 +144 @@
   map<const char *, int, IonTypeCompare> IonTypeLineMap;
   if (LineMapping == NULL) {
-    eLog() << Verbose(0) << "map pointer is NULL: " << LineMapping << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "map pointer is NULL: " << LineMapping << endl);
     performCriticalExit();
     return;
@@ -159 +160 @@
       LineMapping[CurrentLine+(nr++)] = runner->second;
     else {
-      eLog() << Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl);
       performCriticalExit();
     }
@@ -499 +500 @@
 //        case 'j': // BoxLength
 //          Log() << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
+//          double * const cell_size = World::get()->cell_size;
 //          for (int i=0;i<6;i++) {
 //            Log() << Verbose(0) << "Cell size" << i << ": ";
-//            cin >> mol->cell_size[i];
+//            cin >> cell_size[i];
 //          }
 //          break;
@@ -657 +659 @@
 {
   if (FileBuffer != NULL) {
-    eLog() << Verbose(2) << "deleting present FileBuffer in PrepareFileBuffer()." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "deleting present FileBuffer in PrepareFileBuffer()." << endl);
     delete(FileBuffer);
   }
@@ -683 +685 @@
 
   if (mol == NULL) {
-    eLog() << Verbose(0) << "Molecule is not allocated in LoadMolecule(), exit.";
+    DoeLog(0) && (eLog()<< Verbose(0) << "Molecule is not allocated in LoadMolecule(), exit.");
     performCriticalExit();
   }
@@ -689 +691 @@
   ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
   if (MaxTypes == 0) {
-    eLog() << Verbose(0) << "There are no atoms according to MaxTypes in this config file." << endl;
-    performCriticalExit();
+    DoeLog(1) && (eLog()<< Verbose(1) << "There are no atoms according to MaxTypes in this config file." << endl);
+    //performCriticalExit();
   } else {
     // prescan number of ions per type
@@ -708 +710 @@
     sprintf(name,"Ion_Type%i",MaxTypes);
     if (!ParseForParameter(verbose,FileBuffer, (const char*)name, 1, 1, 1, int_type, &value[0], 1, critical)) {
-      eLog() << Verbose(0) << "There are no atoms in the config file!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "There are no atoms in the config file!" << endl);
       performCriticalExit();
       return;
@@ -853 +855 @@
   ifstream *file = new ifstream(filename);
   if (file == NULL) {
-    eLog() << Verbose(1) << "config file " << filename << " missing!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
     return;
   }
@@ -965 +967 @@
   // Unit cell and magnetic field
   ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
-  mol->cell_size[0] = BoxLength[0];
-  mol->cell_size[1] = BoxLength[3];
-  mol->cell_size[2] = BoxLength[4];
-  mol->cell_size[3] = BoxLength[6];
-  mol->cell_size[4] = BoxLength[7];
-  mol->cell_size[5] = BoxLength[8];
+  double * const cell_size = World::get()->cell_size;
+  cell_size[0] = BoxLength[0];
+  cell_size[1] = BoxLength[3];
+  cell_size[2] = BoxLength[4];
+  cell_size[3] = BoxLength[6];
+  cell_size[4] = BoxLength[7];
+  cell_size[5] = BoxLength[8];
   //if (1) fprintf(stderr,"\n");
 
@@ -1062 +1065 @@
       Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
     } else {
-      eLog() << Verbose(1) << "Bond length table loading failed." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
     }
   }
@@ -1091 +1094 @@
   ifstream *file = new ifstream(filename);
   if (file == NULL) {
-    eLog() << Verbose(1) << "config file " << filename << " missing!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
     return;
   }
@@ -1169 +1172 @@
 
   ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
-  mol->cell_size[0] = BoxLength[0];
-  mol->cell_size[1] = BoxLength[3];
-  mol->cell_size[2] = BoxLength[4];
-  mol->cell_size[3] = BoxLength[6];
-  mol->cell_size[4] = BoxLength[7];
-  mol->cell_size[5] = BoxLength[8];
+  double * const cell_size = World::get()->cell_size;
+  cell_size[0] = BoxLength[0];
+  cell_size[1] = BoxLength[3];
+  cell_size[2] = BoxLength[4];
+  cell_size[3] = BoxLength[6];
+  cell_size[4] = BoxLength[7];
+  cell_size[5] = BoxLength[8];
   if (1) fprintf(stderr,"\n");
   config::DoPerturbation = 0;
@@ -1312 +1316 @@
   // bring MaxTypes up to date
   mol->CountElements();
+  const double * const cell_size = World::get()->cell_size;
   ofstream * const output = new ofstream(filename, ios::out);
   if (output != NULL) {
@@ -1382 +1387 @@
     *output << endl;
     *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
-    *output << mol->cell_size[0] << "\t" << endl;
-    *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
-    *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
+    *output << cell_size[0] << "\t" << endl;
+    *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
+    *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
     // FIXME
     *output << endl;
@@ -1428 +1433 @@
     return result;
   } else {
-    eLog() << Verbose(1) << "Cannot open output file:" << filename << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open output file:" << filename << endl);
     return false;
   }
@@ -1450 +1455 @@
     output = new ofstream(fname->str().c_str(), ios::out);
     if (output == NULL) {
-      eLog() << Verbose(1) << "Cannot open mpqc output file:" << fname << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc output file:" << fname << endl);
       delete(fname);
       return false;
@@ -1493 +1498 @@
     output = new ofstream(fname->str().c_str(), ios::out);
     if (output == NULL) {
-      eLog() << Verbose(1) << "Cannot open mpqc hessian output file:" << fname << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc hessian output file:" << fname << endl);
       delete(fname);
       return false;
@@ -1549 +1554 @@
   f = fopen(name, "w" );
   if (f == NULL) {
-    eLog() << Verbose(1) << "Cannot open pdb output file:" << name << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
     return false;
   }
@@ -1604 +1609 @@
   f = fopen(name, "w" );
   if (f == NULL) {
-    eLog() << Verbose(1) << "Cannot open pdb output file:" << name << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
     Free(&elementNo);
     return false;
@@ -1641 +1646 @@
 /** Stores all atoms in a TREMOLO data input file.
  * Note that this format cannot be parsed again.
+ * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
  * \param *filename name of file (without ".in" suffix!)
  * \param *mol pointer to molecule
@@ -1653 +1659 @@
   output = new ofstream(fname->str().c_str(), ios::out);
   if (output == NULL) {
-    eLog() << Verbose(1) << "Cannot open tremolo output file:" << fname << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open tremolo output file:" << fname << endl);
     delete(fname);
     return false;
@@ -1695 +1701 @@
 /** Stores all atoms from all molecules in a TREMOLO data input file.
  * Note that this format cannot be parsed again.
+ * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
  * \param *filename name of file (without ".in" suffix!)
  * \param *MolList pointer to MoleculeListClass containing all atoms
@@ -1707 +1714 @@
   output = new ofstream(fname->str().c_str(), ios::out);
   if (output == NULL) {
-    eLog() << Verbose(1) << "Cannot open tremolo output file:" << fname << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open tremolo output file:" << fname << endl);
     delete(fname);
     return false;
@@ -1747 +1754 @@
       while (Walker->next != (*MolWalker)->end) {
         Walker = Walker->next;
-        *output << AtomNo << "\t";
+        *output << AtomNo+1 << "\t";
         *output << Walker->Name << "\t";
         *output << (*MolWalker)->name << "\t";
-        *output << MolCounter << "\t";
+        *output << MolCounter+1 << "\t";
         *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
         *output << (double)Walker->type->Valence << "\t";
         *output << Walker->type->symbol << "\t";
         for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
-          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ] << "\t";
+          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ]+1 << "\t";
         for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
           *output << "-\t";

src/defs.hpp

@@ -12 +12 @@
 #define MAX_ELEMENTS 128  //!< maximum number of elements for certain lookup tables 
 #define AtomicLengthToAngstroem  0.52917721 //!< conversion factor from atomic length/bohrradius to angstroem
-#define BONDTHRESHOLD 0.5   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii 
 #define AtomicEnergyToKelvin 315774.67  //!< conversion factor from atomic energy to kelvin via boltzmann factor
 #define KelvinToAtomicTemperature 3.1668152e-06    //!< conversion factor for Kelvin to atomic temperature (Hartree over k_B)

src/ellipsoid.cpp

@@ -241 +241 @@
     x = new Vector[PointsToPick];
   } else {
-    eLog() << Verbose(2) << "Given pointer to vector array seems already allocated." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Given pointer to vector array seems already allocated." << endl);
   }
 
@@ -341 +341 @@
     x = new Vector[PointsToPick];
   } else {
-    eLog() << Verbose(2) << "Given pointer to vector array seems already allocated." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Given pointer to vector array seems already allocated." << endl);
   }
 

src/gslvector.cpp

@@ -6 +6 @@
  */
 
+#include <cassert>
+#include <cmath>
+
 #include "gslvector.hpp"
+#include "defs.hpp"
 
 /** Constructor of class GSLVector.
@@ -27 +31 @@
 };
 
+/** Copy constructor of class GSLVector.
+ * Allocates GSL structures and copies components from \a *src.
+ * \param *src source vector
+ */
+GSLVector::GSLVector(const GSLVector & src) : dimension(src.dimension)
+{
+  vector = gsl_vector_alloc(dimension);
+  gsl_vector_memcpy (vector, src.vector);
+};
+
 /** Destructor of class GSLVector.
  * Frees GSL structures
@@ -33 +47 @@
 {
   gsl_vector_free(vector);
-  dimension = 0;
 };
 
@@ -52 +65 @@
  * \return  m-th element of vector
  */
-double GSLVector::Get(size_t m)
+double GSLVector::Get(size_t m) const
 {
   return gsl_vector_get (vector, m);
@@ -72 +85 @@
  * \return pointer to \a m-th element
  */
-double *GSLVector::Pointer(size_t m)
+double *GSLVector::Pointer(size_t m) const
 {
   return gsl_vector_ptr (vector, m);
@@ -82 +95 @@
  * \return const pointer to \a m-th element
  */
-const double *GSLVector::const_Pointer(size_t m)
+const double *GSLVector::const_Pointer(size_t m) const
 {
   return gsl_vector_const_ptr (vector, m);
+};
+
+/** Returns the dimension of the vector.
+ * \return dimension of vector
+ */
+#ifdef HAVE_INLINE
+inline
+#endif
+size_t GSLVector::GetDimension() const
+{
+  return dimension;
 };
 
@@ -128 +152 @@
   return gsl_vector_reverse (vector);
 };
+
+
+/* ========================== Operators =============================== */
+/** Compares GSLVector \a to GSLVector \a b component-wise.
+ * \param a base GSLVector
+ * \param b GSLVector components to add
+ * \return a == b
+ */
+bool operator==(const GSLVector& a, const GSLVector& b)
+{
+  bool status = true;
+  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
+  for (size_t i=0;i<a.GetDimension();i++)
+    status = status && (fabs(a.Get(i) - b.Get(i)) < MYEPSILON);
+  return status;
+};
+
+/** Sums GSLVector \a to this lhs component-wise.
+ * \param a base GSLVector
+ * \param b GSLVector components to add
+ * \return lhs + a
+ */
+const GSLVector& operator+=(GSLVector& a, const GSLVector& b)
+{
+  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
+  for (size_t i=0;i<a.GetDimension();i++)
+    a.Set(i,a.Get(i)+b.Get(i));
+  return a;
+};
+
+/** Subtracts GSLVector \a from this lhs component-wise.
+ * \param a base GSLVector
+ * \param b GSLVector components to add
+ * \return lhs - a
+ */
+const GSLVector& operator-=(GSLVector& a, const GSLVector& b)
+{
+  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
+  for (size_t i=0;i<a.GetDimension();i++)
+    a.Set(i,a.Get(i)-b.Get(i));
+  return a;
+};
+
+/** factor each component of \a a times a double \a m.
+ * \param a base GSLVector
+ * \param m factor
+ * \return lhs.Get(i) * m
+ */
+const GSLVector& operator*=(GSLVector& a, const double m)
+{
+  for (size_t i=0;i<a.GetDimension();i++)
+    a.Set(i,a.Get(i)*m);
+  return a;
+};
+
+/** Sums two GSLVectors \a  and \b component-wise.
+ * \param a first GSLVector
+ * \param b second GSLVector
+ * \return a + b
+ */
+GSLVector const operator+(const GSLVector& a, const GSLVector& b)
+{
+  GSLVector x(a);
+  for (size_t i=0;i<a.GetDimension();i++)
+    x.Set(i,a.Get(i)+b.Get(i));
+  return x;
+};
+
+/** Subtracts GSLVector \a from \b component-wise.
+ * \param a first GSLVector
+ * \param b second GSLVector
+ * \return a - b
+ */
+GSLVector const operator-(const GSLVector& a, const GSLVector& b)
+{
+  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
+  GSLVector x(a);
+  for (size_t i=0;i<a.GetDimension();i++)
+    x.Set(i,a.Get(i)-b.Get(i));
+  return x;
+};
+
+/** Factors given GSLVector \a a times \a m.
+ * \param a GSLVector
+ * \param m factor
+ * \return m * a
+ */
+GSLVector const operator*(const GSLVector& a, const double m)
+{
+  GSLVector x(a);
+  for (size_t i=0;i<a.GetDimension();i++)
+    x.Set(i,a.Get(i)*m);
+  return x;
+};
+
+/** Factors given GSLVector \a a times \a m.
+ * \param m factor
+ * \param a GSLVector
+ * \return m * a
+ */
+GSLVector const operator*(const double m, const GSLVector& a )
+{
+  GSLVector x(a);
+  for (size_t i=0;i<a.GetDimension();i++)
+    x.Set(i,a.Get(i)*m);
+  return x;
+};
+
+ostream& operator<<(ostream& ost, const GSLVector& m)
+{
+  ost << "(";
+  for (size_t i=0;i<m.GetDimension();i++) {
+    ost << m.Get(i);
+    if (i != 2)
+      ost << ",";
+  }
+  ost << ")";
+  return ost;
+};
+
+/* ====================== Checking State ============================ */
+/** Checks whether vector has all components zero.
+ * @return true - vector is zero, false - vector is not
+ */
+bool GSLVector::IsZero() const
+{
+  return (fabs(Get(0))+fabs(Get(1))+fabs(Get(2)) < MYEPSILON);
+};
+
+/** Checks whether vector has length of 1.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool GSLVector::IsOne() const
+{
+  double NormValue = 0.;
+  for (size_t i=dimension;--i;)
+    NormValue += Get(i)*Get(i);
+  return (fabs(NormValue - 1.) < MYEPSILON);
+};
+

src/gslvector.hpp

@@ -18 +18 @@
 #endif
 
+#include <iostream>
 #include <gsl/gsl_vector.h>
 
@@ -32 +33 @@
   GSLVector(size_t m);
   GSLVector(const GSLVector * const src);
+  GSLVector(const GSLVector & src);
   ~GSLVector();
 
   // Accessing
   void SetFromDoubleArray(double *x);
-  double Get(size_t m);
+  double Get(size_t m) const;
   void Set(size_t m, double x);
-  double *Pointer(size_t m);
-  const double *const_Pointer(size_t m);
+  double *Pointer(size_t m) const;
+  const double *const_Pointer(size_t m) const;
+  size_t GetDimension() const;
 
   // Initializing
@@ -50 +53 @@
   int Reverse();
 
+  // checking state
+  bool IsZero() const;
+  bool IsOne() const;
+
 private:
   gsl_vector *vector;
 
-  size_t dimension;
+  const size_t dimension;
 };
+
+ostream & operator << (ostream& ost, const GSLVector &m);
+bool operator==(const GSLVector& a, const GSLVector& b);
+const GSLVector& operator+=(GSLVector& a, const GSLVector& b);
+const GSLVector& operator-=(GSLVector& a, const GSLVector& b);
+const GSLVector& operator*=(GSLVector& a, const double m);
+GSLVector const operator*(const GSLVector& a, const double m);
+GSLVector const operator*(const double m, const GSLVector& a);
+GSLVector const operator+(const GSLVector& a, const GSLVector& b);
+GSLVector const operator-(const GSLVector& a, const GSLVector& b);
+
 
 

src/helpers.hpp

@@ -134 +134 @@
   LookupTable = Calloc<T*>(count, "CreateFatherLookupTable - **LookupTable");
   if (LookupTable == NULL) {
-    eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "LookupTable memory allocation failed!" << endl);
     performCriticalExit();
     status = false;

src/linkedcell.cpp

@@ -46 +46 @@
   min.Zero();
   Log() << Verbose(1) << "Begin of LinkedCell" << endl;
-  if (set->IsEmpty()) {
-    eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
+  if ((set == NULL) || (set->IsEmpty())) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "set is NULL or contains no linked cell nodes!" << endl);
     return;
   }
@@ -79 +79 @@
   Log() << Verbose(2) << "Allocating cells ... ";
   if (LC != NULL) {
-    eLog() << Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
     return;
   }
@@ -122 +122 @@
   Log() << Verbose(1) << "Begin of LinkedCell" << endl;
   if (set->empty()) {
-    eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "set contains no linked cell nodes!" << endl);
     return;
   }
@@ -151 +151 @@
   Log() << Verbose(2) << "Allocating cells ... ";
   if (LC != NULL) {
-    eLog() << Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
     return;
   }
@@ -199 +199 @@
     status = status && ((n[i] >=0) && (n[i] < N[i]));
   if (!status)
-  eLog() << Verbose(1) << "indices are out of bounds!" << endl;
+  DoeLog(1) && (eLog()<< Verbose(1) << "indices are out of bounds!" << endl);
   return status;
 };
@@ -260 +260 @@
     return status;
   } else {
-    eLog() << Verbose(1) << "Node at " << *Walker << " is out of bounds." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Node at " << *Walker << " is out of bounds." << endl);
     return false;
   }

src/log.cpp

@@ -28 +28 @@
 }
 
+/** Checks verbosity for logger.
+ * Is supposed to be used in construct as this:
+ * DoLog(2) && (Log() << Verbose(2) << "message." << endl);
+ * If DoLog does not return true, the right-hand side is not evaluated and we save some time.
+ * \param verbose verbosity level of this message
+ * \return true - print, false - don't
+ */
+bool DoLog(int verbose) {
+  return verbose <= logger::getInstance()->verbosity;
+}
+
+/** Checks verbosity for errorlogger.
+ * Is supposed to be used in construct as this:
+ * DoLog(2) && (Log() << Verbose(2) << "message." << endl);
+ * If DoLog does not return true, the right-hand side is not evaluated and we save some time.
+ * \param verbose verbosity level of this message
+ * \return true - print, false - don't
+ */
+bool DoeLog(int verbose) {
+  return verbose <= errorLogger::getInstance()->verbosity;
+}
+
 /**
  * Prints an error log entry.

src/log.hpp

@@ -15 +15 @@
 class errorLogger * eLog();
 void setVerbosity(int verbosityLevel);
+bool DoLog(int verbose);
+bool DoeLog(int verbose);
 
 #endif /* LOG_HPP_ */

src/molecule.cpp

@@ -23 +23 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "World.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -45 +46 @@
   for(int i=MAX_ELEMENTS;i--;)
     ElementsInMolecule[i] = 0;
-  cell_size[0] = cell_size[2] = cell_size[5]= 20.;
-  cell_size[1] = cell_size[3] = cell_size[4]= 0.;
-  strcpy(name,"none");
+  strcpy(name,World::get()->DefaultName);
 };
 
@@ -159 +158 @@
   double *matrix = NULL;
   bond *Binder = NULL;
+  double * const cell_size = World::get()->cell_size;
 
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -194 +194 @@
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
-    eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
     return false;
     BondRescale = bondlength;
@@ -237 +237 @@
             SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else {
-            eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name;
+            DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name);
           }
         }
@@ -274 +274 @@
       bondangle = TopOrigin->type->HBondAngle[1];
       if (bondangle == -1) {
-        eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+        DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
         return false;
         bondangle = 0;
@@ -396 +396 @@
       break;
     default:
-      eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl);
       AllWentWell = false;
       break;
@@ -447 +447 @@
     Walker->type = elemente->FindElement(shorthand);
     if (Walker->type == NULL) {
-      eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.";
+      DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
       Walker->type = elemente->FindElement(1);
     }
@@ -543 +543 @@
     add(Binder, last);
   } else {
-    eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl);
   }
   return Binder;
@@ -555 +555 @@
 bool molecule::RemoveBond(bond *pointer)
 {
-  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
   pointer->leftatom->RegisterBond(pointer);
   pointer->rightatom->RegisterBond(pointer);
@@ -569 +569 @@
 bool molecule::RemoveBonds(atom *BondPartner)
 {
-  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
   BondList::const_iterator ForeRunner;
   while (!BondPartner->ListOfBonds.empty()) {
@@ -603 +603 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
+  double * const cell_size = World::get()->cell_size;
   cell_size[0] = dim->x[0];
   cell_size[1] = 0.;
@@ -621 +622 @@
     AtomCount--;
   } else
-    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "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--;
@@ -639 +640 @@
     ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
   else
-    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "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--;
@@ -693 +694 @@
 bool molecule::CheckBounds(const Vector *x) const
 {
+  double * const cell_size = World::get()->cell_size;
   bool result = true;
   int j =-1;

src/molecule.hpp

@@ -82 +82 @@
 class molecule : public PointCloud {
   public:
-    double cell_size[6];//!< cell size
     const periodentafel * const elemente; //!< periodic table with each element
     atom *start;        //!< start of atom list

src/molecule_dynamics.cpp

@@ -306 +306 @@
   for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1
     if (Params.DoubleList[i] > 1) {
-      eLog() << Verbose(0) << "Failed to create an injective PermutationMap!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
       performCriticalExit();
     }
@@ -428 +428 @@
             }
             if (Potential > Params.PenaltyConstants[2]) {
-              eLog() << Verbose(1) << "The two-step permutation procedure did not maintain injectivity!" << endl;
+              DoeLog(1) && (eLog()<< Verbose(1) << "The two-step permutation procedure did not maintain injectivity!" << endl);
               exit(255);
             }
             //Log() << Verbose(0) << endl;
           } else {
-            eLog() << Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
+            DoeLog(1) && (eLog()<< Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl);
             exit(255);
           }
@@ -568 +568 @@
     // parse file into ForceMatrix
     if (!Force.ParseMatrix(file, 0,0,0)) {
-      eLog() << Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
       performCriticalExit();
       return false;
     }
     if (Force.RowCounter[0] != AtomCount) {
-      eLog() << Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl);
       performCriticalExit();
       return false;

src/molecule_fragmentation.cpp

@@ -18 +18 @@
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "World.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -112 +113 @@
         testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1)));  // store fragment number and current factor
         if (!testGraphInsert.second) {
-          eLog() << Verbose(0) << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
+          DoeLog(0) && (eLog()<< Verbose(0) << "KeySet file must be corrupt as there are two equal key sets therein!" << endl);
           performCriticalExit();
         }
@@ -213 +214 @@
     Log() << Verbose(0) << "done." << endl;
   } else {
-    eLog() << Verbose(0) << "Unable to open " << line << " for writing keysets!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Unable to open " << line << " for writing keysets!" << endl);
     performCriticalExit();
     status = false;
@@ -371 +372 @@
       }
     } else {
-      eLog() << Verbose(0) << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl);
       performCriticalExit();
     }
@@ -446 +447 @@
     // transmorph graph keyset list into indexed KeySetList
     if (GlobalKeySetList == NULL) {
-      eLog() << Verbose(1) << "Given global key set list (graph) is NULL!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Given global key set list (graph) is NULL!" << endl);
       return false;
     }
@@ -454 +455 @@
     map<int, pair<double,int> > *AdaptiveCriteriaList = ScanAdaptiveFileIntoMap(path, *IndexKeySetList); // (Root No., (Value, Order)) !
     if (AdaptiveCriteriaList->empty()) {
-      eLog() << Verbose(2) << "Unable to parse file, incrementing all." << endl;
+      DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
       while (Walker->next != end) {
         Walker = Walker->next;
@@ -643 +644 @@
         MolecularWalker->Leaf->FragmentBOSSANOVA(FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
       } else {
-        eLog() << Verbose(1) << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
+        DoeLog(1) && (eLog()<< Verbose(1) << "Subgraph " << MolecularWalker << " has no atoms!" << endl);
       }
       FragmentCounter++;  // next fragment list
@@ -847 +848 @@
 
   Leaf->BondDistance = mol->BondDistance;
-  for(int i=NDIM*2;i--;)
-    Leaf->cell_size[i] = mol->cell_size[i];
 
   // first create the minimal set of atoms from the KeySet
@@ -907 +906 @@
       }
     } else {
-      eLog() << Verbose(1) << "Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
     }
     if ((LonelyFlag) && (Leaf->AtomCount > 1)) {
@@ -1142 +1141 @@
         Log() << Verbose(0) << endl;
         //if (!CheckForConnectedSubgraph(FragmentSearch->FragmentSet))
-          //eLog() << Verbose(1) << "The found fragment is not a connected subgraph!" << endl;
+          //DoeLog(1) && (eLog()<< Verbose(1) << "The found fragment is not a connected subgraph!" << endl);
         InsertFragmentIntoGraph(FragmentSearch);
       }
@@ -1658 +1657 @@
   atom *Walker = NULL;
   atom *OtherWalker = NULL;
+  double * const cell_size = World::get()->cell_size;
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   enum Shading *ColorList = NULL;

src/molecule_geometry.cpp

@@ -15 +15 @@
 #include "memoryallocator.hpp"
 #include "molecule.hpp"
+#include "World.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -26 +27 @@
   bool status = true;
   const Vector *Center = DetermineCenterOfAll();
+  double * const cell_size = World::get()->cell_size;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
@@ -46 +48 @@
 {
   bool status = true;
+  double * const cell_size = World::get()->cell_size;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
@@ -226 +229 @@
 void molecule::TranslatePeriodically(const Vector *trans)
 {
+  double * const cell_size = World::get()->cell_size;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
@@ -252 +256 @@
 {
   atom *Walker = start;
+  double * const cell_size = World::get()->cell_size;
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   double *inversematrix = InverseMatrix(cell_size);

src/molecule_graph.cpp

@@ -17 +17 @@
 #include "memoryallocator.hpp"
 #include "molecule.hpp"
+#include "World.hpp"
 
 struct BFSAccounting
@@ -106 +107 @@
   LinkedCell *LC = NULL;
   bool free_BG = false;
+  double * const cell_size = World::get()->cell_size;
 
   if (BG == NULL) {
@@ -493 +495 @@
   bond *Binder = NULL;
 
+  if (AtomCount == 0)
+    return SubGraphs;
   Log() << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
   DepthFirstSearchAnalysis_Init(DFS, this);
@@ -927 +931 @@
     }
     if (i == vertex->ListOfBonds.size()) {
-      eLog() << Verbose(0) << "Error: All Component entries are already occupied!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Error: All Component entries are already occupied!" << endl);
       performCriticalExit();
     }
   } else {
-    eLog() << Verbose(0) << "Error: Given vertex is NULL!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Error: Given vertex is NULL!" << endl);
     performCriticalExit();
   }

src/moleculelist.cpp

@@ -19 +19 @@
 #include "memoryallocator.hpp"
 #include "periodentafel.hpp"
+#include "World.hpp"
 
 /*********************************** Functions for class MoleculeListClass *************************/
@@ -313 +314 @@
   Tesselation *TesselStruct = NULL;
   if ((srcmol == NULL) || (mol == NULL)) {
-    eLog() << Verbose(1) << "Either fixed or variable molecule is given as NULL." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Either fixed or variable molecule is given as NULL." << endl);
     return false;
   }
@@ -321 +322 @@
   FindNonConvexBorder(mol, TesselStruct, (const LinkedCell *&)LCList, 4., NULL);
   if (TesselStruct == NULL) {
-    eLog() << Verbose(1) << "Could not tesselate the fixed molecule." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Could not tesselate the fixed molecule." << endl);
     return false;
   }
@@ -442 +443 @@
     input.open(line.c_str());
     if (input == NULL) {
-      eLog() << Verbose(0) << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << endl << "Unable to open " << line << ", is the directory correct?" << endl);
       performCriticalExit();
       return false;
@@ -638 +639 @@
   int FragmentCounter = 0;
   ofstream output;
-
+  double cell_size_backup[6];
+  double * const cell_size = World::get()->cell_size;
+
+  // backup cell_size
+  for (int i=0;i<6;i++)
+    cell_size_backup[i] = cell_size[i];
   // store the fragments as config and as xyz
   for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
@@ -646 +652 @@
       strcpy(PathBackup, path);
     else {
-      eLog() << Verbose(0) << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl);
       performCriticalExit();
     }
@@ -682 +688 @@
       j += k + 1;
       BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
-      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
+      cell_size[j] = BoxDimension.x[k] * 2.;
     }
     (*ListRunner)->Translate(&BoxDimension);
@@ -724 +730 @@
   // printing final number
   Log() << Verbose(2) << "Final number of fragments: " << FragmentCounter << "." << endl;
+
+  // restore cell_size
+  for (int i=0;i<6;i++)
+    cell_size[i] = cell_size_backup[i];
 
   return result;
@@ -776 +786 @@
 
   // 1. dissect the molecule into connected subgraphs
-  configuration->BG->ConstructBondGraph(mol);
+  if (!configuration->BG->ConstructBondGraph(mol)) {
+    delete (mol);
+    DoeLog(1) && (eLog()<< Verbose(1) << "There are no bonds." << endl);
+    return;
+  }
 
   // 2. scan for connected subgraphs
@@ -783 +797 @@
   Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
   delete(BackEdgeStack);
+  if ((Subgraphs == NULL) || (Subgraphs->next == NULL)) {
+    delete (mol);
+    DoeLog(1) && (eLog()<< Verbose(1) << "There are no atoms." << endl);
+    return;
+  }
 
   // 3. dissect (the following construct is needed to have the atoms not in the order of the DFS, but in
@@ -824 +843 @@
     Walker = mol->start->next;
     if ((Walker->nr <0) || (Walker->nr >= mol->AtomCount)) {
-      eLog() << Verbose(0) << "Index of atom " << *Walker << " is invalid!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Index of atom " << *Walker << " is invalid!" << endl);
       performCriticalExit();
     }
@@ -833 +852 @@
       molecules[FragmentCounter-1]->AddAtom(Walker);    // counting starts at 1
     } else {
-      eLog() << Verbose(0) << "Atom " << *Walker << " not associated to molecule!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Atom " << *Walker << " not associated to molecule!" << endl);
       performCriticalExit();
     }

src/parser.cpp

@@ -160 +160 @@
   //Log() << Verbose(1) << "Opening " << name << " ... "  << input << endl;
   if (input == NULL) {
-    eLog() << Verbose(1) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << endl << "Unable to open " << name << ", is the directory correct?" << endl);
     //performCriticalExit();
     return false;
@@ -183 +183 @@
   //Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
   if (ColumnCounter[MatrixNr] == 0) {
-    eLog() << Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl);
     performCriticalExit();
   }
@@ -199 +199 @@
   //Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
   if (RowCounter[MatrixNr] == 0) {
-    eLog() << Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl);
     performCriticalExit();
   }
@@ -232 +232 @@
     }
   } else {
-    eLog() << Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
+    DoeLog(1) && (eLog()<< Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl); 
   }
   input.close();
@@ -433 +433 @@
               //Log() << Verbose(0) << "Corresponding index in CurrentFragment is " << m << "." << endl;
               if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
-                eLog() << Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
                 performCriticalExit();
                 return false;
@@ -477 +477 @@
     output.open(line.str().c_str(), ios::out);
     if (output == NULL) {
-      eLog() << Verbose(0) << "Unable to open output energy file " << line.str() << "!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Unable to open output energy file " << line.str() << "!" << endl);
       performCriticalExit();
       return false;
@@ -507 +507 @@
   output.open(line.str().c_str(), ios::out);
   if (output == NULL) {
-    eLog() << Verbose(0) << "Unable to open output matrix file " << line.str() << "!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Unable to open output matrix file " << line.str() << "!" << endl);
     performCriticalExit();
     return false;
@@ -664 +664 @@
       int j = Indices[ FragmentNr ][l];
       if (j > RowCounter[MatrixCounter]) {
-        eLog() << Verbose(0) << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
+        DoeLog(0) && (eLog()<< Verbose(0) << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl);
         performCriticalExit();
         return false;
@@ -802 +802 @@
       int j = Indices[ FragmentNr ][l];
       if (j > RowCounter[MatrixCounter]) {
-        eLog() << Verbose(0) << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
+        DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
         performCriticalExit();
         return false;
@@ -810 +810 @@
           int k = Indices[ FragmentNr ][m];
           if (k > ColumnCounter[MatrixCounter]) {
-            eLog() << Verbose(0) << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
+            DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
             performCriticalExit();
             return false;
@@ -863 +863 @@
               //Log() << Verbose(0) << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
               if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
-                eLog() << Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
                 performCriticalExit();
                 return false;
@@ -881 +881 @@
                   //Log() << Verbose(0) << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
                   if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
-                    eLog() << Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                    DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
                     performCriticalExit();
                     return false;

src/periodentafel.cpp

@@ -314 +314 @@
 
   if (!otherstatus)
-    eLog() << Verbose(2) << "Something went wrong while parsing the other databases!" << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Something went wrong while parsing the other databases!" << endl);
 
   return status;

src/stackclass.hpp

@@ -72 +72 @@
     return true;
   } else {
-    eLog() << Verbose(1) << "Stack is full, " << "Stack: CurrentLastEntry " << CurrentLastEntry<< "\tCurrentFirstEntry " << CurrentFirstEntry << "\tNextFreeField " << NextFreeField << "\tEntryCount " << EntryCount << "!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Stack is full, " << "Stack: CurrentLastEntry " << CurrentLastEntry<< "\tCurrentFirstEntry " << CurrentFirstEntry << "\tNextFreeField " << NextFreeField << "\tEntryCount " << EntryCount << "!" << endl);
     return false;
   }
@@ -87 +87 @@
     Walker = StackList[CurrentFirstEntry];
     if (Walker == NULL)
-      eLog() << Verbose(1) << "Stack's field is empty!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Stack's field is empty!" << endl);
     StackList[CurrentFirstEntry] = NULL;
     if (CurrentFirstEntry != CurrentLastEntry) { // hasn't last item been popped as well?
@@ -96 +96 @@
     }
   } else
-    eLog() << Verbose(1) << "Stack is empty!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Stack is empty!" << endl);
   return Walker;
 };
@@ -111 +111 @@
     StackList[CurrentLastEntry] = NULL;
     if (Walker == NULL)
-      eLog() << Verbose(1) << "Stack's field is empty!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Stack's field is empty!" << endl);
     NextFreeField = CurrentLastEntry;
     if (CurrentLastEntry != CurrentFirstEntry)  // has there been more than one item on stack
       CurrentLastEntry = (CurrentLastEntry + (EntryCount-1)) % EntryCount; // step back from current free field to last (modulo does not work in -1, thus go EntryCount-1 instead)
   } else {
-    eLog() << Verbose(1) << "Stack is empty!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Stack is empty!" << endl);
   }
   return Walker;
@@ -151 +151 @@
     } while (i!=NextFreeField);
   else
-    eLog() << Verbose(1) << "Stack is already empty!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Stack is already empty!" << endl);
   if (found) {
     NextFreeField = CurrentLastEntry;

src/tesselation.cpp

@@ -14 +14 @@
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
+#include "triangleintersectionlist.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
@@ -54 +55 @@
   //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
-    eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
   node = NULL;
 };
@@ -187 +188 @@
   }
   if (!triangles.empty())
-    eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
 };
 
@@ -225 +226 @@
   // get the two triangles
   if (triangles.size() != 2) {
-    eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
     return true;
   }
@@ -251 +252 @@
       BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
     else {
-      eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
     }
     node = runner->second->GetThirdEndpoint(this);
@@ -262 +263 @@
       i++;
     } else {
-      eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
       return true;
     }
@@ -367 +368 @@
   }
   if (Counter < 3) {
-    eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl);
     performCriticalExit();
   }
@@ -429 +430 @@
 
   if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
-    eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
     return false;
   }
@@ -474 +475 @@
 };
 
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
- * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+/** Finds the point on the triangle to the point \a *x.
+ * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point.
+ * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the
+ * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down.
  * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
@@ -723 +726 @@
       Runner[i]++;
       if (Runner[i] == endpoints.end()) {
-        eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl;
+        DoeLog(0) && (eLog()<< Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
         performCriticalExit();
       }
@@ -1069 +1072 @@
       runner->second = NULL;
     } else
-      eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);
   }
   Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
@@ -1322 +1325 @@
   else
     {
-      eLog() << Verbose(0) << "No starting triangle found." << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "No starting triangle found." << endl);
     }
 }
@@ -1521 +1524 @@
           TrianglesOnBoundaryCount++;
         } else {
-          eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
+          DoeLog(2) && (eLog()<< Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);
         }
 
@@ -1614 +1617 @@
             if (NewLines[j]->IsConnectedTo(BLS[0])) {
               if (n>2) {
-                eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl;
+                DoeLog(2) && (eLog()<< Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);
                 return false;
               } else
@@ -1631 +1634 @@
       }
     } else { // something is wrong with FindClosestTriangleToPoint!
-      eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "The closest triangle did not produce an intersection!" << endl);
       return false;
     }
@@ -1734 +1737 @@
           OpenLines.erase(CandidateLine);
         } else {
-          eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl;
+          DoeLog(1) && (eLog()<< Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
         }
 
@@ -1840 +1843 @@
       triangle->lines[i] = NULL;  // free'd or not: disconnect
     } else
-      eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "This line " << i << " has already been free'd." << endl);
   }
 
@@ -1894 +1897 @@
       line->endpoints[i] = NULL;  // free'd or not: disconnect
     } else
-      eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);
   }
   if (!line->triangles.empty())
-    eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);
 
   if (LinesOnBoundary.erase(line->Nr))
@@ -2070 +2073 @@
           }
         } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+          DoeLog(1) && (eLog()<< Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
         }
       }
@@ -2214 +2217 @@
 //            if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
 //              // rotated the wrong way!
-//              eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+//              DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
 //            }
 //
@@ -2271 +2274 @@
 //          }
 //        } else {
-//          eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl;
+//          DoeLog(2) && (eLog()<< Verbose(2) << "Baseline is connected to two triangles already?" << endl);
 //        }
 //      } else {
@@ -2278 +2281 @@
 //    }
 //  } else {
-//    eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl;
+//    DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl);
 //  }
 //
@@ -2344 +2347 @@
     if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
-      eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
     }
 
@@ -2355 +2358 @@
 
   if (CandidateLine.pointlist.empty()) {
-    eLog() << Verbose(2) << "Could not find a suitable candidate." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Could not find a suitable candidate." << endl);
     return false;
   }
@@ -2397 +2400 @@
 //          CandidateLine.ShortestAngle = ShortestAngle;
 //        } else {
-////          eLog() << Verbose(1) << "This triangle consisting of ";
+////          DoeLog(1) && (eLog()<< Verbose(1) << "This triangle consisting of ");
 ////          Log() << Verbose(0) << *(*it)->point << ", ";
 ////          Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
@@ -2418 +2421 @@
 //
 //          } else {
-////            eLog() << Verbose(1) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl;
+////            DoeLog(1) && (eLog()<< Verbose(1) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl);
 //            result = false;
 //          }
@@ -2434 +2437 @@
 //    BaseRay = BLS[0];
 //    if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) {
-//      eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl;
+//      DoeLog(1) && (eLog()<< Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl);
 //      exit(255);
 //    }
@@ -2632 +2635 @@
     BaseLineNormal.Zero();
     if (Base->triangles.size() < 2) {
-      eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
       return 0.;
     }
@@ -2671 +2674 @@
   BaseLineNormal.Zero();
   if (Base->triangles.size() < 2) {
-    eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
     return NULL;
   }
@@ -2707 +2710 @@
   // check whether everything is in place to create new lines and triangles
   if (i<4) {
-    eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "We have not gathered enough baselines!" << endl);
     return NULL;
   }
   for (int j=0;j<4;j++)
     if (OldLines[j] == NULL) {
-      eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "We have not gathered enough baselines!" << endl);
       return NULL;
     }
   for (int j=0;j<2;j++)
     if (OldPoints[j] == NULL) {
-      eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "We have not gathered enough endpoints!" << endl);
       return NULL;
     }
@@ -2761 +2764 @@
     Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
   } else {
-    eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);
     return NULL;
   }
@@ -2792 +2795 @@
       N[i] = LC->n[i];
   } else {
-    eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);
     return;
   }
@@ -2932 +2935 @@
     // test whether old center is on the band's plane
     if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl);
       RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
     }
@@ -2942 +2945 @@
       Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
       if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
-        eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+        DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
       }
 
@@ -2951 +2954 @@
         //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
       } else {
-        eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+        DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);
         return;
       }
@@ -2990 +2993 @@
                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
                       if (fabs(radius - otherradius) > HULLEPSILON) {
-                        eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;
+                        DoeLog(1) && (eLog()<< Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl);
                       }
                       // construct both new centers
@@ -3057 +3060 @@
           }
     } else {
-      eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
     }
   } else {
@@ -3116 +3119 @@
 
   if (LinesOnBoundary.empty()) {
-    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl);
     return NULL;
   }
@@ -3143 +3146 @@
           }
         } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+          DoeLog(1) && (eLog()<< Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
         }
       }
@@ -3149 +3152 @@
   // check whether we found some points
   if (points->empty()) {
-    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
     delete(points);
     return NULL;
@@ -3168 +3171 @@
   DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
   if (points == NULL) {
-    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
     return NULL;
   }
@@ -3222 +3225 @@
 };
 
-
 /** Finds the triangle that is closest to a given Vector \a *x.
  * \param *out output stream for debugging
@@ -3235 +3237 @@
         DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
   if (points == NULL) {
-    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
     return NULL;
   }
@@ -3284 +3286 @@
         const double distance = BaseLineIntersection.NormSquared();
         if (Center.NormSquared() > BaseLine.NormSquared()) {
-          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
+          DoeLog(0) && (eLog()<< Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl);
         }
         if ((ClosestLines.empty()) || (distance < MinDistance)) {
@@ -3315 +3317 @@
  * \param *out output stream for debugging
  * \param *x Vector to look from
+ * \param &distance contains found distance on return
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
@@ -3358 +3361 @@
 bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
 {
-  return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON);
-}
+  Info FunctionInfo(__func__);
+  TriangleIntersectionList Intersections(&Point,this,LC);
+
+  return Intersections.IsInside();
+};
 
 /** Returns the distance to the surface given by the tesselation.
@@ -3432 +3438 @@
 };
 
-/** Calculates distance to a tesselated surface.
+/** Calculates minimum distance from \a&Point to a tesselated surface.
  * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
  * \param &Point point to calculate distance from
@@ -3438 +3444 @@
  * \return distance squared to closest point on surface
  */
-double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
-{
-  BoundaryTriangleSet *triangle = FindClosestTriangleToVector(&Point, LC);
-  const double distance = GetDistanceSquaredToTriangle(Point, triangle);
-  return Min(distance, LC->RADIUS);
+double Tesselation::GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  TriangleIntersectionList Intersections(&Point,this,LC);
+
+  return Intersections.GetSmallestDistance();
+};
+
+/** Calculates minimum distance from \a&Point to a tesselated surface.
+ * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
+ * \param &Point point to calculate distance from
+ * \param *LC needed for finding closest points fast
+ * \return distance squared to closest point on surface
+ */
+BoundaryTriangleSet * Tesselation::GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  TriangleIntersectionList Intersections(&Point,this,LC);
+
+  return Intersections.GetClosestTriangle();
 };
 
@@ -3464 +3485 @@
     ReferencePoint = PointRunner->second;
   } else {
-    eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
     ReferencePoint = NULL;
   }
@@ -3496 +3517 @@
 
   if (connectedPoints->empty()) { // if have not found any points
-    eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl);
     return NULL;
   }
@@ -3529 +3550 @@
 
   if (SetOfNeighbours == NULL) {
-    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Could not find any connected points!" << endl);
     delete(connectedCircle);
     return NULL;
@@ -3540 +3561 @@
       PlaneNormal.AddVector(&(*Runner)->NormalVector);
   } else {
-    eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
+    DoeLog(0) && (eLog()<< Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
     performCriticalExit();
   }
@@ -3559 +3580 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
-      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
       performCriticalExit();
     }
@@ -3610 +3631 @@
 
   if (SetOfNeighbours == NULL) {
-    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Could not find any connected points!" << endl);
     delete(connectedCircle);
     return NULL;
@@ -3664 +3685 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
-      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
       performCriticalExit();
     }
@@ -3687 +3708 @@
     InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
     if (!InserterTest.second) {
-      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl);
       performCriticalExit();
     }
@@ -3727 +3748 @@
     ReferencePoint = PointRunner->second;
   } else {
-    eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
     return NULL;
   }
@@ -3744 +3765 @@
       LineRunner = TouchedLine.find(runner->second);
       if (LineRunner == TouchedLine.end()) {
-        eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl;
+        DoeLog(1) && (eLog()<< Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);
       } else if (!LineRunner->second) {
         LineRunner->second = true;
@@ -3774 +3795 @@
                 }
               } else {
-                eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl;
+                DoeLog(1) && (eLog()<< Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl);
                 triangle = NULL;
               }
@@ -3791 +3812 @@
           LineRunner = TouchedLine.find(CurrentLine);
           if (LineRunner == TouchedLine.end())
-            eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl;
+            DoeLog(1) && (eLog()<< Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl);
           else
             LineRunner->second = true;
@@ -3809 +3830 @@
     }
   } else {
-    eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl);
   }
 
@@ -3889 +3910 @@
 
   if (Point == NULL) {
-    eLog() << Verbose(1) << "Point given is NULL." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Point given is NULL." << endl);
   } else {
     // go through its lines and insert all triangles
@@ -3921 +3942 @@
 
   if (point == NULL) {
-    eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);
     return 0.;
   } else
@@ -3931 +3952 @@
   // get list of connected points
   if (point->lines.empty()) {
-    eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);
     return 0.;
   }
@@ -4012 +4033 @@
         MiddleNode = EndNode;
         if (MiddleNode == connectedPath->end()) {
-          eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl;
+          DoeLog(0) && (eLog()<< Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);
           performCriticalExit();
         }
@@ -4031 +4052 @@
         triangle = GetPresentTriangle(TriangleCandidates);
         if (triangle != NULL) {
-          eLog() << Verbose(0) << "New triangle already present, skipping!" << endl;
+          DoeLog(0) && (eLog()<< Verbose(0) << "New triangle already present, skipping!" << endl);
           StartNode++;
           MiddleNode++;
@@ -4069 +4090 @@
           break;
         } else if (connectedPath->size() < 2) { // something's gone wrong!
-          eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl;
+          DoeLog(0) && (eLog()<< Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);
           performCriticalExit();
         } else {
@@ -4222 +4243 @@
     }
     default:
-      eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);
       performCriticalExit();
       break;
@@ -4275 +4296 @@
   // sanity check
   if (LinesOnBoundary.empty()) {
-    eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.";
+    DoeLog(2) && (eLog()<< Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");
     return DegeneratedLines;
   }
@@ -4299 +4320 @@
       Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl;
     else
-      eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);
   }
 
@@ -4459 +4480 @@
     NearestBoundaryPoint = PointRunner->second;
   } else {
-    eLog() << Verbose(1) << "I cannot find the boundary point." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "I cannot find the boundary point." << endl);
     return;
   }
@@ -4527 +4548 @@
     if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
       if (BestLine == BTS->lines[i]){
-        eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl;
+        DoeLog(0) && (eLog()<< Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);
         performCriticalExit();
       }
@@ -4719 +4740 @@
     // connections to either polygon ...
     if (T->size() % 2 != 0) {
-      eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl);
       performCriticalExit();
     }
@@ -4754 +4775 @@
       AddTesselationLine(TPS[1], TPS[2], 2);
       if (TriangleNrs.empty())
-        eLog() << Verbose(0) << "No more free triangle numbers!" << endl;
+        DoeLog(0) && (eLog()<< Verbose(0) << "No more free triangle numbers!" << endl);
       BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...
       AddTesselationTriangle(); // ... and add
@@ -4763 +4784 @@
     }
     if (!TriangleNrs.empty()) {
-      eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl;
+      DoeLog(0) && (eLog()<< Verbose(0) << "There have been less triangles created than removed!" << endl);
     }
     delete(T);  // remove the triangleset

src/tesselation.hpp

@@ -318 +318 @@
     bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
     double GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const;
-    double GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const;
+    double GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const;
+    BoundaryTriangleSet * GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
     DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;

src/tesselationhelpers.cpp

@@ -81 +81 @@
 
   if (fabs(m11) < MYEPSILON)
-    eLog() << Verbose(1) << "three points are colinear." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "three points are colinear." << endl);
 
   center->x[0] =  0.5 * m12/ m11;
@@ -88 +88 @@
 
   if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
-    eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl);
 
   gsl_matrix_free(A);
@@ -192 +192 @@
   //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
   if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
-    eLog() << Verbose(1) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl);
   }
 
@@ -236 +236 @@
   // test whether new center is on the parameter circle's plane
   if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-    eLog() << Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl);
     helper.ProjectOntoPlane(&CirclePlaneNormal);
   }
@@ -242 +242 @@
   // test whether the new center vector has length of CircleRadius
   if (fabs(radius - CircleRadius) > HULLEPSILON)
-    eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
   alpha = helper.Angle(&RelativeOldSphereCenter);
   // make the angle unique by checking the halfplanes/search direction
@@ -512 +512 @@
 //  Vector BaseLineVector, OrthogonalVector, helper;
 //  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-//    eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
+//    DoeLog(1) && (eLog()<< Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl);
 //    //return false;
 //    exit(1);
@@ -764 +764 @@
     }
   } else {
-    eLog() << Verbose(1) << "Given vrmlfile is " << vrmlfile << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Given vrmlfile is " << vrmlfile << "." << endl);
   }
   delete(center);
@@ -839 +839 @@
     *rasterfile << "9\n#  terminating special property\n";
   } else {
-    eLog() << Verbose(1) << "Given rasterfile is " << rasterfile << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Given rasterfile is " << rasterfile << "." << endl);
   }
   IncludeSphereinRaster3D(rasterfile, Tess, cloud);
@@ -951 +951 @@
   // check number of endpoints in *P
   if (P->endpoints.size() != 4) {
-    eLog() << Verbose(1) << "CountTrianglePairContainingPolygon works only on polygons with 4 nodes!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "CountTrianglePairContainingPolygon works only on polygons with 4 nodes!" << endl);
     return 0;
   }
@@ -957 +957 @@
   // check number of triangles in *T
   if (T->size() < 2) {
-    eLog() << Verbose(1) << "Not enough triangles to have pairs!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Not enough triangles to have pairs!" << endl);
     return 0;
   }

src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

@@ -60 +60 @@
 
   // construct molecule (tetraeder of hydrogens) base
+  TestSurfaceMolecule = new molecule(tafel);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(1., 0., 1. );
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(0., 1., 1. );
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(1., 1., 0. );
+  TestSurfaceMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = hydrogen;
+  Walker->node->Init(0., 0., 0. );
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule->AtomCount, 4 );
+
+  TestList = new MoleculeListClass;
+  TestSurfaceMolecule->ActiveFlag = true;
+  TestList->insert(TestSurfaceMolecule);
+
+  // init tesselation and linked cell
+  Surface = new Tesselation;
+  LC = new LinkedCell(TestSurfaceMolecule, 5.);
+  FindNonConvexBorder(TestSurfaceMolecule, Surface, (const LinkedCell *&)LC, 2.5, NULL);
+
+  // add outer atoms
   TestMolecule = new molecule(tafel);
   Walker = new atom();
-  Walker->type = hydrogen;
-  Walker->node->Init(1., 0., 1. );
-  TestMolecule->AddAtom(Walker);
-  Walker = new atom();
-  Walker->type = hydrogen;
-  Walker->node->Init(0., 1., 1. );
-  TestMolecule->AddAtom(Walker);
-  Walker = new atom();
-  Walker->type = hydrogen;
-  Walker->node->Init(1., 1., 0. );
-  TestMolecule->AddAtom(Walker);
-  Walker = new atom();
-  Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
-  TestMolecule->AddAtom(Walker);
-
-  // check that TestMolecule was correctly constructed
-  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
-
-  TestList = new MoleculeListClass;
+  Walker->type = carbon;
+  Walker->node->Init(4., 0., 4. );
+  TestMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = carbon;
+  Walker->node->Init(0., 4., 4. );
+  TestMolecule->AddAtom(Walker);
+  Walker = new atom();
+  Walker->type = carbon;
+  Walker->node->Init(4., 4., 0. );
+  TestMolecule->AddAtom(Walker);
+  // add inner atoms
+  Walker = new atom();
+  Walker->type = carbon;
+  Walker->node->Init(0.5, 0.5, 0.5 );
+  TestMolecule->AddAtom(Walker);
   TestMolecule->ActiveFlag = true;
   TestList->insert(TestMolecule);
-
-  // init tesselation and linked cell
-  Surface = new Tesselation;
-  FindNonConvexBorder(TestMolecule, Surface, (const LinkedCell *&)LC, 2.5, NULL);
-  LC = new LinkedCell(TestMolecule, 5.);
-  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t)6, Surface->LinesOnBoundary.size() );
-  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->TrianglesOnBoundary.size() );
-
-  // add outer atoms
-  Walker = new atom();
-  Walker->type = carbon;
-  Walker->node->Init(4., 0., 4. );
-  TestMolecule->AddAtom(Walker);
-  Walker = new atom();
-  Walker->type = carbon;
-  Walker->node->Init(0., 4., 4. );
-  TestMolecule->AddAtom(Walker);
-  Walker = new atom();
-  Walker->type = carbon;
-  Walker->node->Init(4., 4., 0. );
-  TestMolecule->AddAtom(Walker);
-  // add inner atoms
-  Walker = new atom();
-  Walker->type = carbon;
-  Walker->node->Init(0.5, 0.5, 0.5 );
-  TestMolecule->AddAtom(Walker);
 
   // init maps
@@ -136 +136 @@
 
 
+/** Checks whether setup() does the right thing.
+ */
+void AnalysisCorrelationToSurfaceUnitTest::SurfaceTest()
+{
+  CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule->AtomCount );
+  CPPUNIT_ASSERT_EQUAL( 4, TestMolecule->AtomCount );
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, TestList->ListOfMolecules.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)6, Surface->LinesOnBoundary.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->TrianglesOnBoundary.size() );
+};
+
 void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceTest()
 {
   // do the pair correlation
   surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
   CPPUNIT_ASSERT( surfacemap != NULL );
   CPPUNIT_ASSERT_EQUAL( (size_t)4, surfacemap->size() );
@@ -149 +162 @@
   surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
   // put pair correlation into bins and check with no range
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
   binmap = BinData( surfacemap, 0.5, 0., 0. );
   CPPUNIT_ASSERT_EQUAL( (size_t)1, binmap->size() );
-  //OutputCorrelation ( binmap );
+  OutputCorrelation ( (ofstream *)&cout, binmap );
   tester = binmap->begin();
   CPPUNIT_ASSERT_EQUAL( 0., tester->first );
@@ -162 +176 @@
   BinPairMap::iterator tester;
   surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
   // ... and check with [0., 2.] range
   binmap = BinData( surfacemap, 0.5, 0., 2. );
   CPPUNIT_ASSERT_EQUAL( (size_t)5, binmap->size() );
-  //OutputCorrelation ( binmap );
+//  OutputCorrelation ( (ofstream *)&cout, binmap );
   tester = binmap->begin();
   CPPUNIT_ASSERT_EQUAL( 0., tester->first );
@@ -179 +194 @@
   BinPairMap::iterator tester;
   surfacemap = CorrelationToSurface( TestList, carbon, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
   // put pair correlation into bins and check with no range
   binmap = BinData( surfacemap, 0.5, 0., 0. );
@@ -184 +200 @@
   OutputCorrelation ( (ofstream *)&cout, binmap );
   // inside point is first and must have negative value
-  tester = binmap->lower_bound(2.95); // start depends on the min value and
+  tester = binmap->lower_bound(4.25-0.5); // start depends on the min value and
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 3, tester->second );
   // inner point
-  tester = binmap->lower_bound(-0.5);
+  tester = binmap->lower_bound(0.);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 1, tester->second );
@@ -197 +213 @@
   BinPairMap::iterator tester;
   surfacemap = CorrelationToSurface( TestList, carbon, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
   // ... and check with [0., 2.] range
   binmap = BinData( surfacemap, 0.5, -2., 4. );
   CPPUNIT_ASSERT_EQUAL( (size_t)13, binmap->size() );
-  OutputCorrelation ( (ofstream *)&cout, binmap );
+//  OutputCorrelation ( (ofstream *)&cout, binmap );
   // three outside points
-  tester = binmap->lower_bound(3.);
+  tester = binmap->lower_bound(4.25-0.5);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 3, tester->second );
   // inner point
-  tester = binmap->lower_bound(-0.5);
+  tester = binmap->lower_bound(0.);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 1, tester->second );
-
 };
 

src/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp

@@ -23 +23 @@
 {
     CPPUNIT_TEST_SUITE( AnalysisCorrelationToSurfaceUnitTest ) ;
+    CPPUNIT_TEST ( SurfaceTest );
     CPPUNIT_TEST ( CorrelationToSurfaceTest );
     CPPUNIT_TEST ( CorrelationToSurfaceHydrogenBinNoRangeTest );
@@ -33 +34 @@
       void setUp();
       void tearDown();
+      void SurfaceTest();
       void CorrelationToSurfaceTest();
       void CorrelationToSurfaceHydrogenBinNoRangeTest();
@@ -43 +45 @@
       MoleculeListClass *TestList;
       molecule *TestMolecule;
+      molecule *TestSurfaceMolecule;
       element *hydrogen;
       element *carbon;

src/unittests/gslvectorunittest.cpp

@@ -115 +115 @@
 
 
+/** UnitTest for operators.
+ */
+void GSLVectorTest::OperatorIsTest()
+{
+  GSLVector zero(3);
+  GSLVector unit(3);
+  zero.SetZero();
+  unit.SetZero();
+  unit.Set(1,1.);
+  // summation and scaling
+  CPPUNIT_ASSERT_EQUAL( true, unit.IsOne() );
+  CPPUNIT_ASSERT_EQUAL( false, zero.IsOne() );
+  CPPUNIT_ASSERT_EQUAL( false, unit.IsZero() );
+  CPPUNIT_ASSERT_EQUAL( true, zero.IsZero() );
+};
+
+/** UnitTest for operators.
+ */
+void GSLVectorTest::OperatorAlgebraTest()
+{
+  GSLVector zero(3);
+  GSLVector unit(3);
+  zero.SetZero();
+  unit.SetZero();
+  unit.Set(1,1.);
+  // summation and scaling
+  CPPUNIT_ASSERT_EQUAL( true, (zero+unit).IsOne() );
+  CPPUNIT_ASSERT_EQUAL( true, (zero+unit).IsOne() );
+  CPPUNIT_ASSERT_EQUAL( true, (zero-unit).IsOne() );
+  CPPUNIT_ASSERT_EQUAL( false, (zero-unit).IsZero() );
+  CPPUNIT_ASSERT_EQUAL( true, (zero+zero).IsZero() );
+  CPPUNIT_ASSERT_EQUAL( false, (unit*0.98).IsOne() );
+  CPPUNIT_ASSERT_EQUAL( false, (0.98*unit).IsOne() );
+  CPPUNIT_ASSERT_EQUAL( true, (unit*1.).IsOne() );
+  CPPUNIT_ASSERT_EQUAL( true, (1.*unit).IsOne() );
+
+  CPPUNIT_ASSERT_EQUAL( unit, (zero+unit) );
+  CPPUNIT_ASSERT_EQUAL( zero, (zero+zero) );
+  CPPUNIT_ASSERT_EQUAL( unit, (unit+zero) );
+
+  unit += zero;
+  CPPUNIT_ASSERT_EQUAL( true, unit.IsOne() );
+  unit *= 1.;
+  CPPUNIT_ASSERT_EQUAL( true, unit.IsOne() );
+};
+
 /********************************************** Main routine **************************************/
 

src/unittests/gslvectorunittest.hpp

@@ -22 +22 @@
     CPPUNIT_TEST (CopyTest );
     CPPUNIT_TEST (ExchangeTest );
+    CPPUNIT_TEST (OperatorAlgebraTest );
+    CPPUNIT_TEST (OperatorIsTest );
     CPPUNIT_TEST_SUITE_END();
 
@@ -32 +34 @@
     void CopyTest();
     void ExchangeTest();
+    void OperatorIsTest();
+    void OperatorAlgebraTest();
 
 private:

src/unittests/logunittest.cpp

@@ -43 +43 @@
 
   Log() << Verbose(0) << "Output a log message." << endl;
-  eLog() << Verbose(0) << "Output an error message." << endl;
+  DoeLog(0) && (eLog()<< Verbose(0) << "Output an error message." << endl);
   setVerbosity(3);
   Log() << Verbose(4) << "This should not be printed." << endl;
-  eLog() << Verbose(4) << "This should not be printed." << endl;
+  DoeLog(4) && (eLog()<< Verbose(4) << "This should not be printed." << endl);
 };
 

src/vector.cpp

@@ -15 +15 @@
 #include "vector.hpp"
 #include "verbose.hpp"
+#include "World.hpp"
 
 #include <gsl/gsl_linalg.h>
@@ -26 +27 @@
  */
 Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
+
+/** Constructor of class vector.
+ */
+Vector::Vector(const Vector * const a)
+{
+  x[0] = a->x[0];
+  x[1] = a->x[1];
+  x[2] = a->x[2];
+};
+
+/** Constructor of class vector.
+ */
+Vector::Vector(const Vector &a)
+{
+  x[0] = a.x[0];
+  x[1] = a.x[1];
+  x[2] = a.x[2];
+};
 
 /** Constructor of class vector.
@@ -262 +281 @@
     return true;
   } else {
-    eLog() << Verbose(2) << "Intersection point " << *this << " is not on plane." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Intersection point " << *this << " is not on plane." << endl);
     return false;
   }
 };
 
-/** Calculates the minimum distance of this vector to the plane.
+/** Calculates the minimum distance vector of this vector to the plane.
  * \param *out output stream for debugging
  * \param *PlaneNormal normal of plane
  * \param *PlaneOffset offset of plane
- * \return distance to plane
- */
-double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
+ * \return distance vector onto to plane
+ */
+Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
 {
   Vector temp;
@@ -291 +310 @@
     sign = 0.;
 
-  return (temp.Norm()*sign);
+  temp.Normalize();
+  temp.Scale(sign);
+  return temp;
+};
+
+/** Calculates the minimum distance of this vector to the plane.
+ * \sa Vector::GetDistanceVectorToPlane()
+ * \param *out output stream for debugging
+ * \param *PlaneNormal normal of plane
+ * \param *PlaneOffset offset of plane
+ * \return distance to plane
+ */
+double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
+{
+  return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
 };
 
@@ -783 +816 @@
       x[i] = C.x[i];
   } else {
-    eLog() << Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl);
   }
 };
@@ -835 +868 @@
   x2.SubtractVector(y2);
   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
     return false;
   }
@@ -869 +902 @@
   Zero();
   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
+    DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
     return false;
   }

src/vector.hpp

@@ -27 +27 @@
 
   Vector();
+  Vector(const Vector * const a);
+  Vector(const Vector &a);
   Vector(const double x1, const double x2, const double x3);
   ~Vector();
@@ -67 +69 @@
   void LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors);
   double CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const;
+  Vector GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const;
   bool GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector);
   bool GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *Normal = NULL);