Changeset 2319ed for src

Timestamp:

Aug 3, 2009, 8:10:09 AM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

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

Children:

0dbddc, 357fba

Parents:

e1589e

Message:

We are one step further in fixing the convex hull: There are two functions of Saskia Metzler missing, but then we may proceed with testing whether the simple correction scheme of the convex envelope works, but one thing: Right now we cannot associate a Tesselation to its molecule as the snake bites it's one tail. Hence, the next commit will consist of fixing this bad-OOP issue.

• Makefile.am: Just some alphabetical resorting.
• atom::atom() new copy constructor
• builder.cpp: some output for cluster volume, molecule::AddCopyAtom() uses new copy constructor
• FillBoxWithMolecule() - new function to fill the remainder of the simulation box with some given filler molecules. Makes explicit use of the tesselated surfaces
• find_convex_border() - InsertStraddlingPoints() and CorrectConcaveBaselines() is called to correct for atoms outside the envelope and caused-by concave points
• Tesselation::InsertStraddlingPoints() enlarges the envelope for all atoms found outside, Tesselation::CorrectConcaveBaselines() corrects all found baselines if the adjacent triangles are concave by flipping.
• boundary.cpp: Lots of helper routines for stuff further below:
  • BoundaryLineSet::IsConnectedTo() checks whether lines have common endpoint
  • BoundaryLineSet::CheckConvexityCriterion() checks whether adjacent triangles are convex or concave
  • BoundaryLineSet::ContainsBoundaryPoint() checks whether line contains a given boundary point
  • BoundaryTriangleSet::ContainsBoundaryLine() and BoundaryTriangleSet::ContainsBoundaryPoint() check whether triangle contains boundary line or point
  • BoundaryTriangleSet::GetIntersectionInsideTriangle() calculates intersection point between line and triangle,
• The following routines are needed to check whether point is in- or outside:
  • Vector::GetIntersectionWithPlane() calculates intersection point of a plane and a line.
  • Vector::GetIntersectionOfTwoLinesOnPlane() intersection of two lines on the same plane.
• FIX: Tesselation::AddPoint() - newly created BoundaryPoint is removed if already present.

Problem: We want to associate a Tesselation class with each molecule class. However, so far we have to know about atoms and bond and molecules inside the Tesselation. We have to remove this dependency and create some intermediate class which enables/encapsulates the access to Vectors, e.g. hidden inside the atom class. This is also good OOP! The Tesselation also only needs a set of Vectors, not more!

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• atom.cpp (modified) (1 diff)
• boundary.cpp (modified) (8 diffs)
• boundary.hpp (modified) (4 diffs)
• builder.cpp (modified) (3 diffs)
• molecules.cpp (modified) (1 diff)
• molecules.hpp (modified) (3 diffs)
• vector.cpp (modified) (2 diffs)
• vector.hpp (modified) (1 diff)

src/Makefile.am

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

src/atom.cpp

@@ -30 +30 @@
   MaxOrder = false;
 };
+
+/** Constructor of class atom.
+ *
+ */
+atom::atom(atom *pointer)
+{
+  Name = NULL;
+  previous = NULL;
+  next = NULL;
+  father = this;  // generally, father is itself
+  Ancestor = NULL;
+  type = pointer->type;  // copy element of atom
+  x.CopyVector(&pointer->x); // copy coordination
+  v.CopyVector(&pointer->v); // copy velocity
+  FixedIon = pointer->FixedIon;
+  nr = -1;
+  sort = &nr;
+}
+
 
 /** Destructor of class atom.

src/boundary.cpp

@@ -132 +132 @@
 ;
 
+/** Checks whether we have a common endpoint with given \a *line.
+ * \param *line other line to test
+ * \return true - common endpoint present, false - not connected
+ */
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+{
+  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
+    return true;
+  else
+    return false;
+};
+
+/** Checks whether the adjacent triangles of a baseline are convex or not.
+ * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
+ * If greater/equal M_PI than we are convex.
+ * \param *out output stream for debugging
+ * \return true - triangles are convex, false - concave or less than two triangles connected
+ */
+bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
+{
+  Vector BaseLineNormal;
+  double angle = 0;
+  // get the two triangles
+  if (TrianglesCount != 2) {
+    *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    return false;
+  }
+  // have a normal vector on the base line pointing outwards
+  BaseLineNormal.Zero();
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    BaseLineNormal.AddVector(&runner->second->NormalVector);
+  BaseLineNormal.Normalize();
+  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
+
+  if ((angle - M_PI) > -MYEPSILON)
+    return true;
+  else
+    return false;
+}
+
+/** Checks whether point is any of the two endpoints this line contains.
+ * \param *point point to test
+ * \return true - point is of the line, false - is not
+ */
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<2;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
+
 ostream &
 operator <<(ostream &ost, BoundaryLineSet &a)
@@ -214 +268 @@
 ;
 
-void
-BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+/** Calculates the normal vector for this triangle.
+ * Is made unique by comparison with \a OtherVector to point in the other direction.
+ * \param &OtherVector direction vector to make normal vector unique.
+ */
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
   // get normal vector
@@ -224 +281 @@
   if (NormalVector.Projection(&OtherVector) > 0)
     NormalVector.Scale(-1.);
-}
-;
+};
+
+/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
+ * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
+ * smaller than the first line.
+ * \param *out output stream for debugging
+ * \param *MolCenter offset vector of line
+ * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
+ * \param *Intersection intersection on plane on return
+ * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
+ */
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
+{
+  Vector CrossPoint;
+  Vector helper;
+  int i=0;
+
+  if (Intersection->GetIntersectionWithPlane(out, &NormalVector, &endpoints[0]->node->x, MolCenter, x)) {
+    *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
+    return false;
+  }
+
+  // Calculate cross point between one baseline and the line from the third endpoint to intersection
+  do {
+    CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, &endpoints[i%3]->node->x, &endpoints[(i+1)%3]->node->x, &endpoints[(i+2)%3]->node->x, Intersection);
+    helper.CopyVector(&endpoints[(i+1)%3]->node->x);
+    helper.SubtractVector(&endpoints[i%3]->node->x);
+    i++;
+    if (i>3)
+      break;
+  } while (CrossPoint.NormSquared() < MYEPSILON);
+  if (i>3) {
+    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
+    exit(255);
+  }
+  CrossPoint.SubtractVector(&endpoints[i%3]->node->x);
+
+  // check whether intersection is inside or not by comparing length of intersection and length of cross point
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
+    return true;
+  } else { // outside!
+    Intersection->Zero();
+    return false;
+  }
+};
+
+/** Checks whether lines is any of the three boundary lines this triangle contains.
+ * \param *line line to test
+ * \return true - line is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
+{
+  for(int i=0;i<3;i++)
+    if (line == lines[i])
+      return true;
+  return false;
+};
+
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point point to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<3;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
 
 ostream &
@@ -774 +901 @@
   TesselStruct->TesselateOnBoundary(out, mol);
 
-  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
+  // 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(out, mol))
+    *out << Verbose(1) << "Insertion of straddling points failed!" << endl;
+
+  // 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
+  if (!TesselStruct->CorrectConcaveBaselines(out))
+    *out << Verbose(1) << "Correction of concave baselines failed!" << endl;
+
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -980 +1115 @@
 }
 ;
+
+
+/** Fills the empty space of the simulation box with water/
+ * \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 distance[NDIM] distance between filling molecules in each direction
+ * \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
+ * \return *mol pointer to new molecule with filled atoms
+ */
+molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
+{
+  molecule *Filling = new molecule(filler->elemente);
+  Vector CurrentPosition;
+  int N[NDIM];
+  int n[NDIM];
+  double *M =  filler->ReturnFullMatrixforSymmetric(filler->cell_size);
+  double Rotations[NDIM*NDIM];
+  Vector AtomTranslations;
+  Vector FillerTranslations;
+  Vector FillerDistance;
+  double FillIt = false;
+  atom *Walker = NULL, *Runner = NULL;
+  bond *Binder = NULL;
+
+  // Center filler at origin
+  filler->CenterOrigin(out);
+  filler->Center.Zero();
+
+  // calculate filler grid in [0,1]^3
+  FillerDistance.Init(distance[0], distance[1], distance[2]);
+  FillerDistance.InverseMatrixMultiplication(M);
+  for(int i=0;i<NDIM;i++)
+    N[i] = ceil(1./FillerDistance.x[i]);
+
+  // 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.Init(n[0], n[1], n[2]);
+        CurrentPosition.MatrixMultiplication(M);
+        // Check whether point is in- or outside
+        FillIt = true;
+        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
+          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
+        }
+
+        if (FillIt) {
+          // fill in Filler
+          *out << 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.);
+
+          // go through all atoms
+          Walker = filler->start;
+          while (Walker != filler->end) {
+            Walker = Walker->next;
+            // copy atom ...
+            *Runner = 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) {
+              double phi[NDIM];
+              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)
+              Runner->x.MatrixMultiplication(Rotations);
+            Runner->x.AddVector(&AtomTranslations);
+            Runner->x.AddVector(&FillerTranslations);
+            Runner->x.AddVector(&CurrentPosition);
+            // insert into Filling
+            Filling->AddAtom(Runner);
+          }
+
+          // go through all bonds and add as well
+          Binder = filler->first;
+          while(Binder != filler->last) {
+            Binder = Binder->next;
+          }
+        } else {
+          // leave empty
+          *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
+        }
+      }
+  return Filling;
+};
+
 
 // =========================================================== class TESSELATION ===========================================
@@ -1196 +1443 @@
 
   MolCenter = mol->DetermineCenterOfAll(out);
+  // create a first tesselation with the given BoundaryPoints
   do {
     flag = false;
@@ -1368 +1616 @@
         *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
   } while (flag);
+
+  // exit
   delete(MolCenter);
+};
+
+/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure with atoms
+ * \return true - all straddling points insert, false - something went wrong
+ */
+bool Tesselation::InsertStraddlingPoints(ofstream *out, molecule *mol)
+{
+  Vector Intersection;
+  atom *Walker = mol->start;
+  Vector *MolCenter = mol->DetermineCenterOfAll(out);
+  while (Walker->next != mol->end) {  // we only have to go once through all points, as boundary can become only bigger
+    // get the next triangle
+    BTS = FindClosestTriangleToPoint(out, &Walker->x);
+    if (BTS == NULL) {
+      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
+      return false;
+    }
+    // get the intersection point
+    if (BTS->GetIntersectionInsideTriangle(out, MolCenter, &Walker->x, &Intersection)) {
+      // we have the intersection, check whether in- or outside of boundary
+      if ((MolCenter->DistanceSquared(&Walker->x) - MolCenter->DistanceSquared(&Intersection)) < -MYEPSILON) {
+        // inside, next!
+        *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
+      } else {
+        // outside!
+        *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
+        class BoundaryLineSet *OldLines[3], *NewLines[3];
+        class BoundaryPointSet *OldPoints[3], *NewPoint;
+        // store the three old lines and old points
+        for (int i=0;i<3;i++) {
+          OldLines[i] = BTS->lines[i];
+          OldPoints[i] = BTS->endpoints[i];
+        }
+        // add Walker to boundary points
+        AddPoint(Walker);
+        if (BPS[0] == NULL)
+          NewPoint = BPS[0];
+        else
+          continue;
+        // remove triangle
+        TrianglesOnBoundary.erase(BTS->Nr);
+        // create three new boundary lines
+        for (int i=0;i<3;i++) {
+          BPS[0] = NewPoint;
+          BPS[1] = OldPoints[i];
+          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
+          LinesOnBoundaryCount++;
+        }
+        // create three new triangle with new point
+        for (int i=0;i<3;i++) { // find all baselines
+          BLS[0] = OldLines[i];
+          int n = 1;
+          for (int j=0;j<3;j++) {
+            if (NewLines[j]->IsConnectedTo(BLS[0])) {
+              if (n>2) {
+                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
+                return false;
+              } else
+                BLS[n++] = NewLines[j];
+            }
+          }
+          // create the triangle
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+          TrianglesOnBoundaryCount++;
+        }
+      }
+    } else { // something is wrong with FindClosestTriangleToPoint!
+      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
+      return false;
+    }
+  }
+
+  // exit
+  delete(MolCenter);
+  return true;
+};
+
+/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
+ * \param *out output stream for debugging
+ * \return true - all baselines were corrected, false - there are still concave pieces
+ */
+bool Tesselation::CorrectConcaveBaselines(ofstream *out)
+{
+  class BoundaryTriangleSet *triangle[2];
+  class BoundaryLineSet *OldLines[4], *NewLine;
+  class BoundaryPointSet *OldPoints[2];
+  Vector BaseLineNormal;
+  class BoundaryLineSet *Base = NULL;
+  int OldTriangles[2], OldBaseLine;
+  int i;
+  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
+    Base = baseline->second;
+
+    // check convexity
+    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
+      *out << Verbose(3) << Base << " has two convex triangles." << endl;
+    } else { // not convex!
+      // get the two triangles
+      i=0;
+      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+        triangle[i++] = runner->second;
+      // gather four endpoints and four lines
+      for (int j=0;j<4;j++)
+        OldLines[j] = NULL;
+      for (int j=0;j<2;j++)
+        OldPoints[j] = NULL;
+      i=0;
+      for (int m=0;m<2;m++) { // go over both triangles
+        for (int j=0;j<3;j++) { // all of their endpoints and baselines
+          if (triangle[m]->lines[j] != Base) // pick not the central baseline
+            OldLines[i++] = triangle[m]->lines[j];
+          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
+            OldPoints[m] = triangle[m]->endpoints[j];
+        }
+      }
+      if (i<4) {
+        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+        return false;
+      }
+      for (int j=0;j<4;j++)
+        if (OldLines[j] == NULL) {
+          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+          return false;
+        }
+      for (int j=0;j<2;j++)
+        if (OldPoints[j] == NULL) {
+          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
+          return false;
+        }
+
+      // remove triangles
+      for (int j=0;j<2;j++) {
+        OldTriangles[j] = triangle[j]->Nr;
+        TrianglesOnBoundary.erase(OldTriangles[j]);
+        triangle[j] = NULL;
+      }
+
+      // remove baseline
+      OldBaseLine = Base->Nr;
+      LinesOnBoundary.erase(OldBaseLine);
+      Base = NULL;
+
+      // construct new baseline (with same number as old one)
+      BPS[0] = OldPoints[0];
+      BPS[1] = OldPoints[1];
+      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
+      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
+
+      // construct new triangles with flipped baseline
+      i=-1;
+      if (BLS[0]->IsConnectedTo(OldLines[2]))
+        i=2;
+      if (BLS[0]->IsConnectedTo(OldLines[2]))
+        i=3;
+      if (i!=-1) {
+        BLS[0] = OldLines[0];
+        BLS[1] = OldLines[i];
+        BLS[2] = NewLine;
+        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
+        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
+
+        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
+        BLS[1] = OldLines[1];
+        BLS[2] = NewLine;
+        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
+        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
+      } else {
+        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
+        return false;
+      }
+    }
+  }
+  return true;
+};
+
+
+/** States whether point is in- or outside of a tesselated surface.
+ * \param *pointer point to be checked
+ * \return true - is inside, false - is outside
+ */
+bool Tesselation::IsInside(Vector *pointer)
+{
+
+  // hier kommt dann Saskias Routine hin...
+
+  return true;
+};
+
+/** Finds the closest triangle to a given point.
+ * \param *out output stream for debugging
+ * \param *x second endpoint of line
+ * \return pointer triangle that is closest, NULL if none was found
+ */
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
+{
+  class BoundaryTriangleSet *triangle = NULL;
+
+  // hier kommt dann Saskias Routine hin...
+
+  return triangle;
 };
 
@@ -1382 +1836 @@
   if (InsertUnique.second) // if new point was not present before, increase counter
     PointsOnBoundaryCount++;
+  else {
+    delete(BPS[0]);
+    BPS[0] = NULL;
+  }
 }
 ;

src/boundary.hpp

@@ -54 +54 @@
 
     void AddTriangle(class BoundaryTriangleSet *triangle);
+    bool IsConnectedTo(class BoundaryLineSet *line);
+    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
+    bool CheckConvexityCriterion(ofstream *out);
 
     class BoundaryPointSet *endpoints[2];
@@ -68 +71 @@
 
     void GetNormalVector(Vector &NormalVector);
+    bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
+    bool ContainsBoundaryLine(class BoundaryLineSet *line);
+    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
 
     class BoundaryPointSet *endpoints[3];
@@ -104 +110 @@
     int CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]);
     void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol);
+    class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, Vector *x);
+    bool IsInside(Vector *pointer);
+    bool InsertStraddlingPoints(ofstream *out, molecule *mol);
+    bool CorrectConcaveBaselines(ofstream *out);
 
     PointMap PointsOnBoundary;
@@ -127 +137 @@
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
+molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
 void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename);
 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS);

src/builder.cpp

@@ -54 +54 @@
 #include "helpers.hpp"
 #include "molecules.hpp"
-
 /********************************************* Subsubmenu routine ************************************/
 
@@ -587 +586 @@
         Find_convex_border((ofstream *)&cout, mol, TesselStruct, &LCList, NULL);
         double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, TesselStruct, configuration);
+        cout << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;
         delete(TesselStruct);
       }
@@ -1815 +1815 @@
                 Find_convex_border((ofstream *)&cout, mol, TesselStruct, &LCList, argv[argptr]);
                 double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, TesselStruct, &configuration);
+                cout << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;
                 delete(TesselStruct);
                 argptr+=1;

src/molecules.cpp

@@ -120 +120 @@
 {
   if (pointer != NULL) {
-    atom *walker = new atom();
-    walker->type = pointer->type;  // copy element of atom
-    walker->x.CopyVector(&pointer->x); // copy coordination
-    walker->v.CopyVector(&pointer->v); // copy velocity
-    walker->FixedIon = pointer->FixedIon;
-    walker->sort = &walker->nr;
+    atom *walker = new atom(pointer);
+    walker->Name = (char *) Malloc(sizeof(char)*strlen(pointer->Name)+1, "atom::atom: *Name");
+    strcpy (walker->Name, pointer->Name);
     walker->nr = last_atom++;  // increase number within molecule
-    walker->father = pointer; //->GetTrueFather();
-    walker->Name = (char *) Malloc(sizeof(char)*strlen(pointer->Name)+1, "molecule::AddCopyAtom: *Name");
-    strcpy (walker->Name, pointer->Name);
     add(walker, end);
     if ((pointer->type != NULL) && (pointer->type->Z != 1))

src/molecules.hpp

@@ -35 +35 @@
 class config;
 class molecule;
+class MoleculeLeafClass;
 class MoleculeListClass;
 class Verbose;
+class Tesselation;
 
 /******************************** Some definitions for easier reading **********************************/
@@ -148 +150 @@
 
   atom();
+  atom(class atom *pointer);
   ~atom();
 
@@ -199 +202 @@
 
 ostream & operator << (ostream &ost, const bond &b);
-
-class MoleculeLeafClass;
-
 
 #define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 

src/vector.cpp

@@ -196 +196 @@
 };
 
+/** Calculates the intersection point between a line defined by \a *LineVector and \a *LineVector2 and a plane defined by \a *Normal and \a *PlaneOffset.
+ * According to [Bronstein] the vectorial plane equation is:
+ *   -# \f$\stackrel{r}{\rightarrow} \cdot \stackrel{N}{\rightarrow} + D = 0\f$,
+ * where \f$\stackrel{r}{\rightarrow}\f$ is the vector to be testet, \f$\stackrel{N}{\rightarrow}\f$ is the plane's normal vector and
+ * \f$D = - \stackrel{a}{\rightarrow} \stackrel{N}{\rightarrow}\f$, the offset with respect to origin, if \f$\stackrel{a}{\rightarrow}\f$,
+ * is an offset vector onto the plane. The line is parametrized by \f$\stackrel{x}{\rightarrow} + k \stackrel{t}{\rightarrow}\f$, where
+ * \f$\stackrel{x}{\rightarrow}\f$ is the offset and \f$\stackrel{t}{\rightarrow}\f$ the directional vector (NOTE: No need to normalize
+ * the latter). Inserting the parametrized form into the plane equation and solving for \f$k\f$, which we insert then into the parametrization
+ * of the line yields the intersection point on the plane.
+ * \param *out output stream for debugging
+ * \param *PlaneNormal Plane's normal vector
+ * \param *PlaneOffset Plane's offset vector
+ * \param *LineVector first vector of line
+ * \param *LineVector2 second vector of line
+ * \return true -  \a this contains intersection point on return, false - line is parallel to plane
+ */
+bool Vector::GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *LineVector, Vector *LineVector2)
+{
+  double factor;
+  Vector Direction;
+
+  // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
+  Direction.CopyVector(LineVector2);
+  Direction.SubtractVector(LineVector);
+  if (Direction.ScalarProduct(PlaneNormal) < MYEPSILON)
+    return false;
+  factor = LineVector->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
+  Direction.Scale(factor);
+  CopyVector(LineVector);
+  SubtractVector(&Direction);
+
+  // test whether resulting vector really is on plane
+  Direction.CopyVector(this);
+  Direction.SubtractVector(PlaneOffset);
+  if (Direction.ScalarProduct(PlaneNormal) > MYEPSILON)
+    return true;
+  else
+    return false;
+};
+
+/** Calculates the intersection of the two lines that are both on the same plane.
+ * Note that we do not check whether they are on the same plane.
+ * \param *out output stream for debugging
+ * \param *Line1a first vector of first line
+ * \param *Line1b second vector of first line
+ * \param *Line2a first vector of second line
+ * \param *Line2b second vector of second line
+ * \return true - \a this will contain the intersection on return, false - lines are parallel
+ */
+bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b)
+{
+  double k1,a1,h1,b1,k2,a2,h2,b2;
+  // equation for line 1
+  if (fabs(Line1a->x[0] - Line2a->x[0]) < MYEPSILON) {
+    k1 = 0;
+    h1 = 0;
+  } else {
+    k1 = (Line1a->x[1] - Line2a->x[1])/(Line1a->x[0] - Line2a->x[0]);
+    h1 = (Line1a->x[2] - Line2a->x[2])/(Line1a->x[0] - Line2a->x[0]);
+  }
+  a1 = 0.5*((Line1a->x[1] + Line2a->x[1]) - k1*(Line1a->x[0] + Line2a->x[0]));
+  b1 = 0.5*((Line1a->x[2] + Line2a->x[2]) - h1*(Line1a->x[0] + Line2a->x[0]));
+
+  // equation for line 2
+  if (fabs(Line2a->x[0] - Line2a->x[0]) < MYEPSILON) {
+    k1 = 0;
+    h1 = 0;
+  } else {
+    k1 = (Line2a->x[1] - Line2a->x[1])/(Line2a->x[0] - Line2a->x[0]);
+    h1 = (Line2a->x[2] - Line2a->x[2])/(Line2a->x[0] - Line2a->x[0]);
+  }
+  a1 = 0.5*((Line2a->x[1] + Line2a->x[1]) - k1*(Line2a->x[0] + Line2a->x[0]));
+  b1 = 0.5*((Line2a->x[2] + Line2a->x[2]) - h1*(Line2a->x[0] + Line2a->x[0]));
+
+  // calculate cross point
+  if (fabs((a1-a2)*(h1-h2) - (b1-b2)*(k1-k2)) < MYEPSILON) {
+    x[0] = (a2-a1)/(k1-k2);
+    x[1] = (k1*a2-k2*a1)/(k1-k2);
+    x[2] = (h1*b2-h2*b1)/(h1-h2);
+    *out << Verbose(4) << "Lines do intersect at " << this << "." << endl;
+    return true;
+  } else {
+    *out << Verbose(4) << "Lines do not intersect." << endl;
+    return false;
+  }
+};
+
 /** Calculates the projection of a vector onto another \a *y.
  * \param *y array to second vector
@@ -453 +540 @@
 };
 
-/** Do a matrix multiplication with \a *matrix' inverse.
+/** Do a matrix multiplication with the \a *A' inverse.
  * \param *matrix NDIM_NDIM array
  */

src/vector.hpp

@@ -46 +46 @@
   void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors);
   double CutsPlaneAt(Vector *A, Vector *B, Vector *C);
+  bool GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *LineVector, Vector *LineVector2);
+  bool GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b);
   bool GetOneNormalVector(const Vector *x1);
   bool MakeNormalVector(const Vector *y1);