source: src/Tesselation/boundary.cpp@ d103d3

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2011 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/** \file boundary.cpp
 *
 * Implementations and super-function for envelopes
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "atom.hpp"
#include "Bond/bond.hpp"
#include "boundary.hpp"
#include "BoundaryLineSet.hpp"
#include "BoundaryPointSet.hpp"
#include "BoundaryTriangleSet.hpp"
#include "Box.hpp"
#include "CandidateForTesselation.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "linkedcell.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "PointCloudAdaptor.hpp"
#include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
#include "RandomNumbers/RandomNumberGenerator.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "World.hpp"

#include <iostream>
#include <iomanip>

#include<gsl/gsl_poly.h>
#include<time.h>

// ========================================== F U N C T I O N S =================================


/** Determines greatest diameters of a cluster defined by its convex envelope.
 * Looks at lines parallel to one axis and where they intersect on the projected planes
 * \param *out output stream for debugging
 * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
 * \param *mol molecule structure representing the cluster
 * \param *&TesselStruct Tesselation structure with triangles
 * \param IsAngstroem whether we have angstroem or atomic units
 * \return NDIM array of the diameters
 */
double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
{
        Info FunctionInfo(__func__);
  // get points on boundary of NULL was given as parameter
  bool BoundaryFreeFlag = false;
  double OldComponent = 0.;
  double tmp = 0.;
  double w1 = 0.;
  double w2 = 0.;
  Vector DistanceVector;
  Vector OtherVector;
  int component = 0;
  int Othercomponent = 0;
  Boundaries::const_iterator Neighbour;
  Boundaries::const_iterator OtherNeighbour;
  double *GreatestDiameter = new double[NDIM];

  const Boundaries *BoundaryPoints;
  if (BoundaryPtr == NULL) {
    BoundaryFreeFlag = true;
    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  } else {
    BoundaryPoints = BoundaryPtr;
    LOG(0, "Using given boundary points set.");
  }
  // determine biggest "diameter" of cluster for each axis
  for (int i = 0; i < NDIM; i++)
    GreatestDiameter[i] = 0.;
  for (int axis = 0; axis < NDIM; axis++)
    { // regard each projected plane
      //LOG(1, "Current axis is " << axis << ".");
      for (int j = 0; j < 2; j++)
        { // and for both axis on the current plane
          component = (axis + j + 1) % NDIM;
          Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
          //LOG(1, "Current component is " << component << ", Othercomponent is " << Othercomponent << ".");
          for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
              //LOG(1, "Current runner is " << *(runner->second.second) << ".");
              // seek for the neighbours pair where the Othercomponent sign flips
              Neighbour = runner;
              Neighbour++;
              if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                Neighbour = BoundaryPoints[axis].begin();
              DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
              do { // seek for neighbour pair where it flips
                  OldComponent = DistanceVector[Othercomponent];
                  Neighbour++;
                  if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                    Neighbour = BoundaryPoints[axis].begin();
                  DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
                  //LOG(2, "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << ".");
                } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                  OldComponent) - DistanceVector[Othercomponent] / fabs(
                  DistanceVector[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
              if (runner != Neighbour) {
                  OtherNeighbour = Neighbour;
                  if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
                    OtherNeighbour = BoundaryPoints[axis].end();
                  OtherNeighbour--;
                  //LOG(1, "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << ".");
                  // now we have found the pair: Neighbour and OtherNeighbour
                  OtherVector = (runner->second.second->getPosition()) - (OtherNeighbour->second.second->getPosition());
                  //LOG(1, "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << ".");
                  //LOG(1, "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << ".");
                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                  w1 = fabs(OtherVector[Othercomponent]);
                  w2 = fabs(DistanceVector[Othercomponent]);
                  tmp = fabs((w1 * DistanceVector[component] + w2
                      * OtherVector[component]) / (w1 + w2));
                  // mark if it has greater diameter
                  //LOG(1, "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << ".");
                  GreatestDiameter[component] = (GreatestDiameter[component]
                      > tmp) ? GreatestDiameter[component] : tmp;
                } //else
              //LOG(1, "Saw no sign flip, probably top or bottom node.");
            }
        }
    }
  LOG(0, "RESULT: The biggest diameters are "
      << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
      << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
      : "atomiclength") << ".");

  // free reference lists
  if (BoundaryFreeFlag)
    delete[] (BoundaryPoints);

  return GreatestDiameter;
}
;


/** Determines the boundary points of a cluster.
 * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
 * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
 * center and first and last point in the triple, it is thrown out.
 * \param *out output stream for debugging
 * \param *mol molecule structure representing the cluster
 * \param *&TesselStruct pointer to Tesselation structure
 */
Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
{
        Info FunctionInfo(__func__);
  PointMap PointsOnBoundary;
  LineMap LinesOnBoundary;
  TriangleMap TrianglesOnBoundary;
  Vector *MolCenter = mol->DetermineCenterOfAll();
  Vector helper;
  BoundariesTestPair BoundaryTestPair;
  Vector AxisVector;
  Vector AngleReferenceVector;
  Vector AngleReferenceNormalVector;
  Vector ProjectedVector;
  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, Nr)
  double angle = 0.;

  // 3a. Go through every axis
  for (int axis = 0; axis < NDIM; axis++) {
    AxisVector.Zero();
    AngleReferenceVector.Zero();
    AngleReferenceNormalVector.Zero();
    AxisVector[axis] = 1.;
    AngleReferenceVector[(axis + 1) % NDIM] = 1.;
    AngleReferenceNormalVector[(axis + 2) % NDIM] = 1.;

    LOG(1, "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << ".");

    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
      ProjectedVector = (*iter)->getPosition() - (*MolCenter);
      ProjectedVector.ProjectOntoPlane(AxisVector);

      // correct for negative side
      const double radius = ProjectedVector.NormSquared();
      if (fabs(radius) > MYEPSILON)
        angle = ProjectedVector.Angle(AngleReferenceVector);
      else
        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues

      //LOG(1, "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << ".");
      if (ProjectedVector.ScalarProduct(AngleReferenceNormalVector) > 0) {
        angle = 2. * M_PI - angle;
      }
    LOG(1, "Inserting " << **iter << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector);
      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, TesselPointDistancePair (radius, (*iter))));
      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
        LOG(2, "Encountered two vectors whose projection onto axis " << axis << " is equal: ");
        LOG(2, "Present vector: " << *BoundaryTestPair.first->second.second);
        LOG(2, "New vector: " << **iter);
        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
          BoundaryTestPair.first->second.second = (*iter);
          LOG(2, "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << ".");
        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
          helper = (*iter)->getPosition() - (*MolCenter);
          const double oldhelperNorm = helper.NormSquared();
          helper = BoundaryTestPair.first->second.second->getPosition() - (*MolCenter);
          if (helper.NormSquared() < oldhelperNorm) {
            BoundaryTestPair.first->second.second = (*iter);
            LOG(2, "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << ".");
          } else {
            LOG(2, "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << ".");
          }
        } else {
          LOG(2, "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << ".");
        }
      }
    }
    // printing all inserted for debugging
    //    {
    //      std::stringstream output;
    //      output << "Printing list of candidates for axis " << axis << " which we have inserted so far: ";
    //      int i=0;
    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    //        if (runner != BoundaryPoints[axis].begin())
    //          output << ", " << i << ": " << *runner->second.second;
    //        else
    //          output << i << ": " << *runner->second.second;
    //        i++;
    //      }
    //      LOG(1, output.str());
    //    }
    // 3c. throw out points whose distance is less than the mean of left and right neighbours
    bool flag = false;
    LOG(1, "Looking for candidates to kick out by convex condition ... ");
    do { // do as long as we still throw one out per round
      flag = false;
      Boundaries::iterator left = BoundaryPoints[axis].begin();
      Boundaries::iterator right = BoundaryPoints[axis].begin();
      Boundaries::iterator runner = BoundaryPoints[axis].begin();
      bool LoopOnceDone = false;
      while (!LoopOnceDone) {
        runner = right;
        right++;
        // set neighbours correctly
        if (runner == BoundaryPoints[axis].begin()) {
          left = BoundaryPoints[axis].end();
        } else {
          left = runner;
        }
        left--;
        if (right == BoundaryPoints[axis].end()) {
          right = BoundaryPoints[axis].begin();
          LoopOnceDone = true;
        }
        // check distance

        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
        {
          Vector SideA, SideB, SideC, SideH;
          SideA = left->second.second->getPosition() - (*MolCenter);
          SideA.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideA: " << SideA);

          SideB = right->second.second->getPosition() -(*MolCenter);
          SideB.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideB: " << SideB);

          SideC = left->second.second->getPosition() - right->second.second->getPosition();
          SideC.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideC: " << SideC);

          SideH = runner->second.second->getPosition() -(*MolCenter);
          SideH.ProjectOntoPlane(AxisVector);
          //          LOG(1, "SideH: " << SideH);

          // calculate each length
          const double a = SideA.Norm();
          //const double b = SideB.Norm();
          //const double c = SideC.Norm();
          const double h = SideH.Norm();
          // calculate the angles
          const double alpha = SideA.Angle(SideH);
          const double beta = SideA.Angle(SideC);
          const double gamma = SideB.Angle(SideH);
          const double delta = SideC.Angle(SideH);
          const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
          //LOG(1, " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << ".");
          LOG(1, "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << ".");
          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
            // throw out point
            LOG(1, "Throwing out " << *runner->second.second << ".");
            BoundaryPoints[axis].erase(runner);
            runner = right;
            flag = true;
          }
        }
      }
    } while (flag);
  }
  delete(MolCenter);
  return BoundaryPoints;
};

/** Tesselates the convex boundary by finding all boundary points.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's.
 * \param *BoundaryPts set of boundary points to use or NULL
 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList atoms in LinkedCell list
 * \param *filename filename prefix for output of vertex data
 * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
 */
void FindConvexBorder(const molecule* mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
{
        Info FunctionInfo(__func__);
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

  if (TesselStruct != NULL) // free if allocated
    delete(TesselStruct);
  TesselStruct = new class Tesselation;

  // 1. Find all points on the boundary
  if (BoundaryPts == NULL) {
    BoundaryFreeFlag = true;
    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  } else {
    BoundaryPoints = BoundaryPts;
    LOG(0, "Using given boundary points set.");
  }

// printing all inserted for debugging
  if (DoLog(1)) {
    for (int axis=0; axis < NDIM; axis++) {
      std::stringstream output;
      output << "Printing list of candidates for axis " << axis << " which we have inserted so far: ";
      int i=0;
      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
        if (runner != BoundaryPoints[axis].begin())
          output << ", " << i << ": " << *runner->second.second;
        else
          output << i << ": " << *runner->second.second;
        i++;
      }
      LOG(1, output.str());
    }
  }

  // 2. fill the boundary point list
  for (int axis = 0; axis < NDIM; axis++)
    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
          LOG(2, "Point " << *(runner->second.second) << " is already present.");

  LOG(0, "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary.");
  // now we have the whole set of edge points in the BoundaryList

  // listing for debugging
  //if (DoLog(1)) {
  //  std::stringstream output;
  //  output << "Listing PointsOnBoundary:";
  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
  //    output << " " << *runner->second;
  //  }
  //  LOG(1, output.str());
  //}

  // 3a. guess starting triangle
  TesselStruct->GuessStartingTriangle();

  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
  PointCloudAdaptor< molecule > cloud(const_cast<molecule *>(mol), mol->name);
  TesselStruct->TesselateOnBoundary(cloud);

  // 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(cloud, LCList))
    ELOG(1, "Insertion of straddling points failed!");

  LOG(0, "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points.");

  // 4. Store triangles in tecplot file
  StoreTrianglesinFile(mol, TesselStruct, filename, "_intermed");

  // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
  bool AllConvex = true;
  class BoundaryLineSet *line = NULL;
  do {
    AllConvex = true;
    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
      line = LineRunner->second;
      LOG(1, "INFO: Current line is " << *line << ".");
      if (!line->CheckConvexityCriterion()) {
        LOG(1, "... line " << *line << " is concave, flipping it.");

        // flip the line
        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
          ELOG(1, "Correction of concave baselines failed!");
        else {
          TesselStruct->FlipBaseline(line);
          LOG(1, "INFO: Correction of concave baselines worked.");
          LineRunner = TesselStruct->LinesOnBoundary.begin(); // LineRunner may have been erase if line was deleted from LinesOnBoundary
        }
      }
    }
  } while (!AllConvex);

  // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
//  if (!TesselStruct->CorrectConcaveTesselPoints(out))
//    ELOG(1, "Correction of concave tesselpoints failed!");

  LOG(0, "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points.");

  // 4. Store triangles in tecplot file
  StoreTrianglesinFile(mol, TesselStruct, filename, "");

  // free reference lists
  if (BoundaryFreeFlag)
    delete[] (BoundaryPoints);
};

/** For testing removes one boundary point after another to check for leaks.
 * \param *out output stream for debugging
 * \param *TesselStruct Tesselation containing envelope with boundary points
 * \param *mol molecule
 * \param *filename name of file
 * \return true - all removed, false - something went wrong
 */
bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
{
        Info FunctionInfo(__func__);
  int i=0;
  char number[MAXSTRINGSIZE];

  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    ELOG(1, "TesselStruct is empty.");
    return false;
  }

  PointMap::iterator PointRunner;
  while (!TesselStruct->PointsOnBoundary.empty()) {
    if (DoLog(1)) {
      std::stringstream output;
      output << "Remaining points are: ";
      for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
        output << *(PointSprinter->second) << "\t";
      LOG(1, output.str());
    }

    PointRunner = TesselStruct->PointsOnBoundary.begin();
    // remove point
    TesselStruct->RemovePointFromTesselatedSurface(PointRunner->second);

    // store envelope
    sprintf(number, "-%04d", i++);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, number);
  }

  return true;
};

/** Creates a convex envelope from a given non-convex one.
 * -# First step, remove concave spots, i.e. singular "dents"
 *   -# We go through all PointsOnBoundary.
 *   -# We CheckConvexityCriterion() for all its lines.
 *   -# If all its lines are concave, it cannot be on the convex envelope.
 *   -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
 *   -# We calculate the additional volume.
 *   -# We go over all lines until none yields a concavity anymore.
 * -# Second step, remove concave lines, i.e. line-shape "dents"
 *   -# We go through all LinesOnBoundary
 *   -# We CheckConvexityCriterion()
 *   -# If it returns concave, we flip the line in this quadrupel of points (abusing the degeneracy of the tesselation)
 *   -# We CheckConvexityCriterion(),
 *   -# if it's concave, we continue
 *   -# if not, we mark an error and stop
 * Note: This routine - for free - calculates the difference in volume between convex and
 * non-convex envelope, as the former is easy to calculate - VolumeOfConvexEnvelope() - it
 * can be used to compute volumes of arbitrary shapes.
 * \param *out output stream for debugging
 * \param *TesselStruct non-convex envelope, is changed in return!
 * \param *mol molecule
 * \param *filename name of file
 * \return volume difference between the non- and the created convex envelope
 */
double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
{
        Info FunctionInfo(__func__);
  double volume = 0;
  class BoundaryPointSet *point = NULL;
  class BoundaryLineSet *line = NULL;
  bool Concavity = false;
  char dummy[MAXSTRINGSIZE];
  PointMap::iterator PointRunner;
  PointMap::iterator PointAdvance;
  LineMap::iterator LineRunner;
  LineMap::iterator LineAdvance;
  TriangleMap::iterator TriangleRunner;
  TriangleMap::iterator TriangleAdvance;
  int run = 0;

  // check whether there is something to work on
  if (TesselStruct == NULL) {
    ELOG(1, "TesselStruct is empty!");
    return volume;
  }

  // First step: RemovePointFromTesselatedSurface
  do {
    Concavity = false;
    sprintf(dummy, "-first-%d", run);
    //CalculateConcavityPerBoundaryPoint(TesselStruct);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);

    PointRunner = TesselStruct->PointsOnBoundary.begin();
    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
    while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
      PointAdvance++;
      point = PointRunner->second;
      LOG(1, "INFO: Current point is " << *point << ".");
      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
        line = LineRunner->second;
        LOG(1, "INFO: Current line of point " << *point << " is " << *line << ".");
        if (!line->CheckConvexityCriterion()) {
          // remove the point if needed
          LOG(1, "... point " << *point << " cannot be on convex envelope.");
          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
          sprintf(dummy, "-first-%d", ++run);
          StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
          Concavity = true;
          break;
        }
      }
      PointRunner = PointAdvance;
    }

    sprintf(dummy, "-second-%d", run);
    //CalculateConcavityPerBoundaryPoint(TesselStruct);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);

    // second step: PickFarthestofTwoBaselines
    LineRunner = TesselStruct->LinesOnBoundary.begin();
    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
      LineAdvance++;
      line = LineRunner->second;
      LOG(1, "INFO: Picking farthest baseline for line is " << *line << ".");
      // take highest of both lines
      if (TesselStruct->IsConvexRectangle(line) == NULL) {
        const double tmp = TesselStruct->PickFarthestofTwoBaselines(line);
        volume += tmp;
        if (tmp != 0.) {
          TesselStruct->FlipBaseline(line);
          Concavity = true;
        }
      }
      LineRunner = LineAdvance;
    }
    run++;
  } while (Concavity);
  //CalculateConcavityPerBoundaryPoint(TesselStruct);
  //StoreTrianglesinFile(mol, filename, "-third");

  // third step: IsConvexRectangle
//  LineRunner = TesselStruct->LinesOnBoundary.begin();
//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
//    LineAdvance++;
//    line = LineRunner->second;
//    LOG(1, "INFO: Current line is " << *line << ".");
//    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
//      //LOG(1, "INFO: Next line will be " << *(LineAdvance->second) << ".");
//    if (!line->CheckConvexityCriterion(out)) {
//      LOG(1, "INFO: ... line " << *line << " is concave, flipping it.");
//
//      // take highest of both lines
//      point = TesselStruct->IsConvexRectangle(line);
//      if (point != NULL)
//        volume += TesselStruct->RemovePointFromTesselatedSurface(point);
//    }
//    LineRunner = LineAdvance;
//  }

  CalculateConcavityPerBoundaryPoint(TesselStruct);
  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");

  // end
  LOG(0, "Volume is " << volume << ".");
  return volume;
};


/** Determines the volume of a cluster.
 * Determines first the convex envelope, then tesselates it and calculates its volume.
 * \param *out output stream for debugging
 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *configuration needed for path to store convex envelope file
 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
 */
double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration)
{
        Info FunctionInfo(__func__);
  bool IsAngstroem = configuration->GetIsAngstroem();
  double volume = 0.;
  Vector x;
  Vector y;

  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
    { // go through every triangle, calculate volume of its pyramid with CoG as peak
      x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1);
      y = runner->second->getEndpoint(0) - runner->second->getEndpoint(2);
      const double a = x.Norm();
      const double b = y.Norm();
      const double c = runner->second->getEndpoint(2).distance(runner->second->getEndpoint(1));
      const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
      x = runner->second->getPlane().getNormal();
      x.Scale(runner->second->getEndpoint(1).ScalarProduct(x));
      const double h = x.Norm(); // distance of CoG to triangle
      const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
      LOG(1, "INFO: Area of triangle is " << setprecision(10) << G << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
          << h << " and the volume is " << PyramidVolume << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
      volume += PyramidVolume;
    }
  LOG(0, "RESULT: The summed volume is " << setprecision(6)
      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");

  return volume;
};

/** Stores triangles to file.
 * \param *out output stream for debugging
 * \param *mol molecule with atoms and bonds
 * \param *TesselStruct Tesselation with boundary triangles
 * \param *filename prefix of filename
 * \param *extraSuffix intermediate suffix
 */
void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
{
        Info FunctionInfo(__func__);
  PointCloudAdaptor< molecule > cloud(const_cast<molecule *>(mol), mol->name);
  // 4. Store triangles in tecplot file
  if (filename != NULL) {
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      WriteTecplotFile(tecplot, TesselStruct, cloud, -1);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      WriteRaster3dFile(rasterplot, TesselStruct, cloud);
      rasterplot->close();
      delete(rasterplot);
    }
  }
};

/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
 * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
 * TODO: Here, we need a VolumeOfGeneralEnvelope (i.e. non-convex one)
 * \param *out output stream for debugging
 * \param *configuration needed for path to store convex envelope file
 * \param *mol molecule structure representing the cluster
 * \param *&TesselStruct Tesselation structure with triangles on return
 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
 * \param celldensity desired average density in final cell
 */
void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
{
        Info FunctionInfo(__func__);
  bool IsAngstroem = true;
  double *GreatestDiameter = NULL;
  Boundaries *BoundaryPoints = NULL;
  class Tesselation *TesselStruct = NULL;
  Vector BoxLengths;
  int repetition[NDIM] = { 1, 1, 1 };
  int TotalNoClusters = 1;
  double totalmass = 0.;
  double clustervolume = 0.;
  double cellvolume = 0.;

  // transform to PAS by Action
  Vector MainAxis(0.,0.,1.);
  mol->RotateToPrincipalAxisSystem(MainAxis);

  IsAngstroem = configuration->GetIsAngstroem();
  BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
  PointCloudAdaptor< molecule > cloud(mol, mol->name);
  LinkedCell *LCList = new LinkedCell(cloud, 10.);
  FindConvexBorder(mol, BoundaryPoints, TesselStruct, (const LinkedCell *&)LCList, NULL);
  delete (LCList);
  delete[] BoundaryPoints;


  // some preparations beforehand
  if (ClusterVolume == 0)
    clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
  else
    clustervolume = ClusterVolume;

  delete TesselStruct;

  for (int i = 0; i < NDIM; i++)
    TotalNoClusters *= repetition[i];

  // sum up the atomic masses
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
      totalmass += (*iter)->getType()->getMass();
  }
  LOG(0, "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit.");
  LOG(0, "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");

  // solve cubic polynomial
  LOG(1, "Solving equidistant suspension in water problem ...");
  if (IsAngstroem)
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
  else
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
  LOG(1, "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");

  double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
  LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
  if (minimumvolume > cellvolume) {
    ELOG(1, "the containing box already has a greater volume than the envisaged cell volume!");
    LOG(0, "Setting Box dimensions to minimum possible, the greatest diameters.");
    for (int i = 0; i < NDIM; i++)
      BoxLengths[i] = GreatestDiameter[i];
    mol->CenterEdge(&BoxLengths);
  } else {
    BoxLengths[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
    BoxLengths[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
    BoxLengths[2] = minimumvolume - cellvolume;
    double x0 = 0.;
    double x1 = 0.;
    double x2 = 0.;
    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
      LOG(0, "RESULT: The resulting spacing is: " << x0 << " .");
    else {
      LOG(0, "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " .");
      x0 = x2; // sorted in ascending order
    }

    cellvolume = 1.;
    for (int i = 0; i < NDIM; i++) {
      BoxLengths[i] = repetition[i] * (x0 + GreatestDiameter[i]);
      cellvolume *= BoxLengths[i];
    }

    // set new box dimensions
    LOG(0, "Translating to box with these boundaries.");
    mol->SetBoxDimension(&BoxLengths);
    mol->CenterInBox();
  }
  delete GreatestDiameter;
  // update Box of atoms by boundary
  mol->SetBoxDimension(&BoxLengths);
  LOG(0, "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
};


/** Fills the empty space around other molecules' surface of the simulation box with a filler.
 * \param *out output stream for debugging
 * \param *List list of molecules already present in box
 * \param *TesselStruct contains tesselated surface
 * \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
 * \param epsilon distance to surface which is not filled
 * \param RandAtomDisplacement maximum distance for random displacement per atom
 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
 * \param DoRandomRotation true - do random rotiations, false - don't
 */
void 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__);
  molecule *Filling = World::getInstance().createMolecule();
  Vector CurrentPosition;
  int N[NDIM];
  int n[NDIM];
  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
  RealSpaceMatrix Rotations;
  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
  Vector AtomTranslations;
  Vector FillerTranslations;
  Vector FillerDistance;
  Vector Inserter;
  double FillIt = false;
  bond *Binder = NULL;
  double phi[NDIM];
  map<molecule *, Tesselation *> TesselStruct;
  map<molecule *, LinkedCell *> LCList;

  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
    if ((*ListRunner)->getAtomCount() > 0) {
      LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
      PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
      LCList[(*ListRunner)] = new LinkedCell(cloud, 10.); // get linked cell list
      LOG(1, "Pre-creating tesselation for molecule " << *ListRunner << ".");
      TesselStruct[(*ListRunner)] = NULL;
      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
    }

  // Center filler at origin
  filler->CenterEdge(&Inserter);
  const int FillerCount = filler->getAtomCount();
  LOG(2, "INFO: Filler molecule has the following bonds:");
  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner) {
      if ((*BondRunner)->leftatom == *AtomRunner)
        LOG(2, "  " << *(*BondRunner));
    }
  }

  atom * CopyAtoms[FillerCount];

  // calculate filler grid in [0,1]^3
  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
  for(int i=0;i<NDIM;i++)
    N[i] = (int) ceil(1./FillerDistance[i]);
  LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");

  // initialize seed of random number generator to current time
  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
  const double rng_min = random.min();
  const double rng_max = random.max();
  //srand ( time(NULL) );

  // go over [0,1]^3 filler grid
  for (n[0] = 0; n[0] < N[0]; n[0]++)
    for (n[1] = 0; n[1] < N[1]; n[1]++)
      for (n[2] = 0; n[2] < N[2]; n[2]++) {
        // calculate position of current grid vector in untransformed box
        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
        // create molecule random translation vector ...
        for (int i=0;i<NDIM;i++)
          FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
        LOG(2, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");

        // go through all atoms
        for (int i=0;i<FillerCount;i++)
          CopyAtoms[i] = NULL;

        // have same rotation angles for all molecule's atoms
        if (DoRandomRotation)
          for (int i=0;i<NDIM;i++)
            phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);

        for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){

          // create atomic random translation vector ...
          for (int i=0;i<NDIM;i++)
            AtomTranslations[i] = RandomAtomDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);

          // ... and rotation matrix
          if (DoRandomRotation) {
            Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
            Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
            Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
            Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
            Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
            Rotations.set(1,2,              sin(phi[1])                                                    );
            Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
            Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
            Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
          }

          // ... and put at new position
          Inserter = (*iter)->getPosition();
          if (DoRandomRotation)
            Inserter *= Rotations;
          Inserter += AtomTranslations + FillerTranslations + CurrentPosition;

          // check whether inserter is inside box
          Inserter *= MInverse;
          FillIt = true;
          for (int i=0;i<NDIM;i++)
            FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
          Inserter *= 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(1, "INFO: Position at " << Inserter << " is outer point.");
            // copy atom ...
            CopyAtoms[(*iter)->getNr()] = (*iter)->clone();
            (*CopyAtoms[(*iter)->getNr()]).setPosition(Inserter);
            Filling->AddAtom(CopyAtoms[(*iter)->getNr()]);
            LOG(1, "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->getNr()]->getPosition()) << ".");
          } else {
            LOG(1, "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance.");
            CopyAtoms[(*iter)->getNr()] = NULL;
            continue;
          }
        }
        // go through all bonds and add as well
        for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
          const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
          for(BondList::const_iterator BondRunner = ListOfBonds.begin();
              BondRunner != ListOfBonds.end();
              ++BondRunner)
            if ((*BondRunner)->leftatom == *AtomRunner) {
              Binder = (*BondRunner);
              if ((CopyAtoms[Binder->leftatom->getNr()] != NULL) && (CopyAtoms[Binder->rightatom->getNr()] != NULL)) {
                LOG(3, "Adding Bond between " << *CopyAtoms[Binder->leftatom->getNr()] << " and " << *CopyAtoms[Binder->rightatom->getNr()]<< ".");
                Filling->AddBond(CopyAtoms[Binder->leftatom->getNr()], CopyAtoms[Binder->rightatom->getNr()], Binder->BondDegree);
              }
            }
        }
      }
  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    delete LCList[*ListRunner];
    delete TesselStruct[(*ListRunner)];
  }
};

/** Rotates given molecule \a Filling and moves its atoms according to given
 *  \a RandomAtomDisplacement.
 *
 *  Note that for rotation to be sensible, the molecule should be centered at
 *  the origin. This is not done here!
 *
 *  \param &Filling molecule whose atoms to displace
 *  \param RandomAtomDisplacement magnitude of random displacement
 *  \param &Rotations 3D rotation matrix (or unity if no rotation desired)
 */
void RandomizeMoleculePositions(
    molecule *&Filling,
    double RandomAtomDisplacement,
    RealSpaceMatrix &Rotations,
    RandomNumberGenerator &random
    )
{
  const double rng_min = random.min();
  const double rng_max = random.max();

  Vector AtomTranslations;
  for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
    Vector temp = (*miter)->getPosition();
    temp *= Rotations;
    (*miter)->setPosition(temp);
    // create atomic random translation vector ...
    for (int i=0;i<NDIM;i++)
      AtomTranslations[i] = RandomAtomDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
    (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
  }
}

/** Removes all atoms of a molecule outside.
 *
 * If the molecule is empty, it is removed as well.
 *
 * @param Filling molecule whose atoms to check, removed if eventually left
 *        empty.
 * @return true - atoms had to be removed, false - no atoms have been removed
 */
bool RemoveAtomsOutsideDomain(molecule *&Filling)
{
  bool status = false;
  Box &Domain = World::getInstance().getDomain();
  // check if all is still inside domain
  for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
    // check whether each atom is inside box
    if (!Domain.isInside((*miter)->getPosition())) {
      status = true;
      atom *Walker = *miter;
      ++miter;
      World::getInstance().destroyAtom(Walker);
    } else {
      ++miter;
    }
  }
  if (Filling->empty()) {
    LOG(0, "Removing molecule " << Filling->getName() << ", all atoms have been removed.");
    World::getInstance().destroyMolecule(Filling);
    Filling = NULL;
  }
  return status;
}

/** Checks whether there are no atoms inside a sphere around \a CurrentPosition
 *  except those atoms present in \a *filler.
 *  If filler is NULL, then we just call LinkedCell::GetPointsInsideSphere() and
 *  check whether the return list is empty.
 * @param *filler
 * @param boundary
 * @param CurrentPosition
 */
bool isSpaceAroundPointVoid(
    LinkedCell *LC,
    molecule *filler,
    const double boundary,
    Vector &CurrentPosition)
{
  size_t compareTo = 0;
  TesselPointSTLList* liste = LC->GetPointsInsideSphere(boundary == 0. ? MYEPSILON : boundary, &CurrentPosition);
  if (filler != NULL) {
    for (TesselPointSTLList::const_iterator iter = liste->begin();
        iter != liste->end();
        ++iter) {
      for (molecule::iterator miter = filler->begin();
          miter != filler->end();
          ++miter) {
        if (*iter == *miter)
          ++compareTo;
      }
    }
  }
  const bool result = (liste->size() == compareTo);
  if (!result) {
    LOG(0, "Skipping because of the following atoms:");
    for (TesselPointSTLList::const_iterator iter = liste->begin();
        iter != liste->end();
        ++iter) {
      LOG(0, **iter);
    }
  }
  delete(liste);
  return result;
}

/** Sets given 3x3 matrix to a random rotation matrix.
 *
 * @param a matrix to set
 */
inline void setRandomRotation(RealSpaceMatrix &a)
{
  double phi[NDIM];
  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
  const double rng_min = random.min();
  const double rng_max = random.max();

  for (int i=0;i<NDIM;i++) {
    phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);
    std::cout << "Random angle is " << phi[i] << std::endl;
  }

  a.setRotation(phi);
}

/** Fills the empty space of the simulation box with water.
 * \param *filler molecule which the box is to be filled with
 * \param configuration contains box dimensions
 * \param distance[NDIM] distance between filling molecules in each direction
 * \param boundary length of boundary zone between molecule and filling molecules
 * \param RandAtomDisplacement maximum distance for random displacement per atom
 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
 * \param MinDistance minimum distance to boundary of domain and present molecules
 * \param DoRandomRotation true - do random rotations, false - don't
 */
void FillVoidWithMolecule(
    molecule *&filler,
    config &configuration,
    const double distance[NDIM],
    const double boundary,
    const double RandomAtomDisplacement,
    const double RandomMolDisplacement,
    const double MinDistance,
    const bool DoRandomRotation
    )
{
  Info FunctionInfo(__func__);
  molecule *Filling = NULL;
  Vector CurrentPosition;
  int N[NDIM];
  int n[NDIM];
  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
  RealSpaceMatrix Rotations;
  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
  Vector FillerTranslations;
  Vector FillerDistance;
  Vector Inserter;
  double FillIt = false;
  Vector firstInserter;
  bool firstInsertion = true;
  const Box &Domain = World::getInstance().getDomain();
  map<molecule *, LinkedCell *> LCList;
  std::vector<molecule *> List = World::getInstance().getAllMolecules();
  MoleculeListClass *MolList = World::getInstance().getMolecules();

  for (std::vector<molecule *>::iterator ListRunner = List.begin(); ListRunner != List.end(); ListRunner++)
    if ((*ListRunner)->getAtomCount() > 0) {
      LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
      PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
      LCList[(*ListRunner)] = new LinkedCell(cloud, 10.); // get linked cell list
    }

  // Center filler at its center of gravity
  Vector *gravity = filler->DetermineCenterOfGravity();
  filler->CenterAtVector(gravity);
  delete gravity;
  //const int FillerCount = filler->getAtomCount();
  LOG(2, "INFO: Filler molecule has the following bonds:");
  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner)
        LOG(2, "  " << *(*BondRunner));
  }

  // calculate filler grid in [0,1]^3
  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
  for(int i=0;i<NDIM;i++)
    N[i] = (int) ceil(1./FillerDistance[i]);
  LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");

  // initialize seed of random number generator to current time
  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
  const double rng_min = random.min();
  const double rng_max = random.max();
  //srand ( time(NULL) );

  // go over [0,1]^3 filler grid
  for (n[0] = 0; n[0] < N[0]; n[0]++)
    for (n[1] = 0; n[1] < N[1]; n[1]++)
      for (n[2] = 0; n[2] < N[2]; n[2]++) {
        // calculate position of current grid vector in untransformed box
        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
        // create molecule random translation vector ...
        for (int i=0;i<NDIM;i++) // have the random values [-1,1]*RandomMolDisplacement
          FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
        LOG(2, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");

        // ... and rotation matrix
        if (DoRandomRotation)
          setRandomRotation(Rotations);
        else
          Rotations.setIdentity();


        // Check whether there is anything too close by and whether atom is outside of domain
        FillIt = true;
        for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
          FillIt = FillIt && isSpaceAroundPointVoid(
              ListRunner->second,
              (firstInsertion ? filler : NULL),
              boundary,
              CurrentPosition);
          FillIt = FillIt && (Domain.isInside(CurrentPosition))
              && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
          if (!FillIt)
            break;
        }

        // insert into Filling
        if (FillIt) {
          Inserter = CurrentPosition + FillerTranslations;
          LOG(1, "INFO: Position at " << Inserter << " is void point.");
          // fill!
          Filling = filler->CopyMolecule();
          RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations, random);
          // translation
          Filling->Translate(&Inserter);
          // remove out-of-bounds atoms
          const bool status = RemoveAtomsOutsideDomain(Filling);
          if ((firstInsertion) && (!status)) { // no atom has been removed
            // remove copied atoms and molecule again
            Filling->removeAtomsinMolecule();
            World::getInstance().destroyMolecule(Filling);
            // and mark is final filler position
            Filling = filler;
            firstInsertion = false;
            firstInserter = Inserter;
          } else {
            // TODO: Remove when World has no MoleculeListClass anymore
            if (Filling)
              MolList->insert(Filling);
          }
        } else {
         LOG(1, "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance.");
          continue;
        }
      }

  // have we inserted any molecules?
  if (firstInsertion) {
    // If not remove filler
    for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
      atom *Walker = *miter;
      World::getInstance().destroyAtom(Walker);
    }
    World::getInstance().destroyMolecule(filler);
  } else {
    // otherwise translate and randomize to final position
    if (DoRandomRotation)
      setRandomRotation(Rotations);
    else
      Rotations.setIdentity();
    RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations, random);
    // translation
    filler->Translate(&firstInserter);
    // remove out-of-bounds atoms
    RemoveAtomsOutsideDomain(filler);
  }

  LOG(0, MolList->ListOfMolecules.size() << " molecules have been inserted.");

  for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
    delete ListRunner->second;
    LCList.erase(ListRunner);
  }
};

/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's
 * \param *&TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *&LCList atoms in LinkedCell list
 * \param RADIUS radius of the virtual sphere
 * \param *filename filename prefix for output of vertex data
 * \return true - tesselation successful, false - tesselation failed
 */
bool FindNonConvexBorder(molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
{
        Info FunctionInfo(__func__);
  bool freeLC = false;
  bool status = false;
  CandidateForTesselation *baseline = NULL;
  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
  bool TesselationFailFlag = false;

  mol->getAtomCount();

  if (TesselStruct == NULL) {
    LOG(1, "Allocating Tesselation struct ...");
    TesselStruct= new Tesselation;
  } else {
    delete(TesselStruct);
    LOG(1, "Re-Allocating Tesselation struct ...");
    TesselStruct = new Tesselation;
  }

  // initialise Linked Cell
  PointCloudAdaptor< molecule > cloud(mol, mol->name);
  if (LCList == NULL) {
    LCList = new LinkedCell(cloud, 2.*RADIUS);
    freeLC = true;
  }

  // 1. get starting triangle
  if (!TesselStruct->FindStartingTriangle(RADIUS, LCList)) {
    ELOG(0, "No valid starting triangle found.");
    //performCriticalExit();
  }
  if (filename != NULL) {
    if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
      TesselStruct->Output(filename, cloud);
    }
  }

  // 2. expand from there
  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
    (cerr << "There are " <<  TesselStruct->TrianglesOnBoundary.size() << " triangles and " << TesselStruct->OpenLines.size() << " open lines to scan for candidates." << endl);
    // 2a. print OpenLines without candidates
    LOG(1, "There are the following open lines to scan for a candidates:");
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
      if (Runner->second->pointlist.empty())
        LOG(1, " " << *(Runner->second));

    // 2b. find best candidate for each OpenLine
    TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);

    // 2c. print OpenLines with candidates again
    LOG(1, "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:");
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
      LOG(1, " " << *(Runner->second));

    // 2d. search for smallest ShortestAngle among all candidates
    double ShortestAngle = 4.*M_PI;
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
      if (Runner->second->ShortestAngle < ShortestAngle) {
        baseline = Runner->second;
        ShortestAngle = baseline->ShortestAngle;
        LOG(1, "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle);
      }
    }
    // 2e. if we found one, add candidate
    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
      OneLoopWithoutSuccessFlag = false;
    else {
      TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
    }

    // 2f. write temporary envelope
    if (filename != NULL) {
      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
        TesselStruct->Output(filename, cloud);
      }
    }
  }
//  // check envelope for consistency
//  status = CheckListOfBaselines(TesselStruct);
//
//  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
//  //->InsertStraddlingPoints(mol, LCList);
//  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
//  class TesselPoint *Runner = NULL;
//    Runner = *iter;
//    LOG(1, "Checking on " << Runner->Name << " ... ");
//    if (!->IsInnerPoint(Runner, LCList)) {
//      LOG(2, Runner->Name << " is outside of envelope, adding via degenerated triangles.");
//      ->AddBoundaryPointByDegeneratedTriangle(Runner, LCList);
//    } else {
//      LOG(2, Runner->Name << " is inside of or on envelope.");
//    }
//  }

//  // Purges surplus triangles.
//  TesselStruct->RemoveDegeneratedTriangles();
//
//  // check envelope for consistency
//  status = CheckListOfBaselines(TesselStruct);

  cout << "before correction" << endl;

  // store before correction
  StoreTrianglesinFile(mol, TesselStruct, filename, "");

//  // correct degenerated polygons
//  TesselStruct->CorrectAllDegeneratedPolygons();
//
  // check envelope for consistency
  status = CheckListOfBaselines(TesselStruct);

  // write final envelope
  CalculateConcavityPerBoundaryPoint(TesselStruct);
  cout << "after correction" << endl;
  StoreTrianglesinFile(mol, TesselStruct, filename, "");

  if (freeLC)
    delete(LCList);

  return status;
};


/** Finds a hole of sufficient size in \a *mols to embed \a *srcmol into it.
 * \param *out output stream for debugging
 * \param *mols molecules in the domain to embed in between
 * \param *srcmol embedding molecule
 * \return *Vector new center of \a *srcmol for embedding relative to \a this
 */
Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
{
        Info FunctionInfo(__func__);
  Vector *Center = new Vector;
  Center->Zero();
  // calculate volume/shape of \a *srcmol

  // find embedding holes

  // if more than one, let user choose

  // return embedding center
  return Center;
};