Changeset 49de64 for src

Timestamp:

Jun 10, 2008, 11:21:32 AM (17 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:

cdee6b

Parents:

233e33

Message:

VolumeOfConvexEnvelope() finished and tested with SiOH4, 1_2_dimethoxyethane and CSHCluster (Ratio1-1)

Location:

src

Files:

• molecules.cpp (modified) (9 diffs)
• molecules.hpp (modified) (2 diffs)

src/molecules.cpp

@@ -1381 +1381 @@
         
         // 1. calculate center of gravity
+        *out << endl;
         vector *CenterOfGravity = DetermineCenterOfGravity(out);
         
+        // 2. translate all points into CoG
+        *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+        Walker=start;
+        while (Walker->next != end) {
+          Walker = Walker->next;
+          Walker->x.Translate(CenterOfGravity);
+        }
         // 2. calculate distance of each atom and sort into hash table
 //      map<double, int> DistanceSet;
@@ -1393 +1401 @@
         
         // 3. Find all points on the boundary
+  *out << Verbose(1) << "Finding all boundary points." << endl;
         Boundaries BoundaryPoints[NDIM];        // first is alpha, second is (r, nr)
         BoundariesTestPair BoundaryTestPair;
-        vector AxisVector;
+        vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+        double radius, angle;
         // 3a. Go through every axis
         for (int axis=0; axis<NDIM; axis++)  {
                 AxisVector.Zero();
+    AngleReferenceVector.Zero();
+    AngleReferenceNormalVector.Zero();
                 AxisVector.x[axis] = 1.;
+                AngleReferenceVector.x[(axis+1)%NDIM] = 1.;
+                AngleReferenceNormalVector.x[(axis+2)%NDIM] = 1.;
+//              *out << Verbose(1) << "Axisvector is ";
+//              AxisVector.Output(out);
+//              *out << " and AngleReferenceVector is ";
+//              AngleReferenceVector.Output(out);
+//              *out << "." << endl;
+//    *out << " and AngleReferenceNormalVector is ";
+//    AngleReferenceNormalVector.Output(out);
+//    *out << "." << endl;
                 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
                 Walker = start;
@@ -1406 +1428 @@
                         vector ProjectedVector;
                         ProjectedVector.CopyVector(&Walker->x);
-                        ProjectedVector.Scale(-1.*Walker->x.Projection(&AxisVector));
-                        ProjectedVector.AddVector(&Walker->x);  // subtract projection from vector
-                        BoundaryTestPair = BoundaryPoints[axis].insert( BoundariesPair (ProjectedVector.Angle(&AxisVector), DistanceNrPair (ProjectedVector.Norm(), Walker) ) );
+                        ProjectedVector.ProjectOntoPlane(&AxisVector);
+                  // correct for negative side
+                  //if (Projection(y) < 0)
+                    //angle = 2.*M_PI - angle;
+                        radius = ProjectedVector.Norm();
+                        if (fabs(radius) > MYEPSILON)
+                          angle = ProjectedVector.Angle(&AngleReferenceVector);
+                        else
+                          angle = 0.;  // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+                          
+                        //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+                        if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+                          angle = 2.*M_PI - angle;
+                        }
+                        *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
+                        ProjectedVector.Output(out);
+                        *out << endl;
+                        BoundaryTestPair = BoundaryPoints[axis].insert( BoundariesPair (angle, DistanceNrPair (radius, Walker) ) );
                         if (BoundaryTestPair.second) { // successfully inserted
-                        } else { // same point exists, check first r, then distance of original vectors to centers of gravity
+                        } else { // same point exists, check first r, then distance of original vectors to center of gravity
                                 *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
                                 *out << Verbose(2) << "Present vector: ";
@@ -1424 +1461 @@
                                         *out << Verbose(2) << "Keeping new vector." << endl;
                                 } else if (tmp == BoundaryTestPair.first->second.first) {
-                                        if (BoundaryTestPair.first->second.second->x.Distance(CenterOfGravity) < Walker->x.Distance(CenterOfGravity)) {
+                                        if (BoundaryTestPair.first->second.second->x.ScalarProduct(&BoundaryTestPair.first->second.second->x) < Walker->x.ScalarProduct(&Walker->x)) { // Norm() does a sqrt, which makes it a lot slower
                                                 BoundaryTestPair.first->second.second = Walker;
                                                 *out << Verbose(2) << "Keeping new vector." << endl;
@@ -1435 +1472 @@
                         }
                 }
-                // 3c. throw out points whose distance is less than the mean of left and right neighbours
+                // printing all inserted for debugging
+                {
+      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      int i=0;
+      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        if (runner != BoundaryPoints[axis].begin())
+          *out << ", " << i << ": " << *runner->second.second;
+        else
+          *out << i << ": " << *runner->second.second;
+        i++;
+      }
+      *out << endl;
+                }
+    // 3c. throw out points whose distance is less than the mean of left and right neighbours
                 bool flag = false;
                 do { // do as long as we still throw one out per round
+                  *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
                         flag = false;
                         Boundaries::iterator left = BoundaryPoints[axis].end();
@@ -1455 +1506 @@
                                 }
                                 // check distance
-                                if (runner->second.first < (left->second.first + right->second.first)/2. ) {
-                                        // throw out point
-                                        BoundaryPoints[axis].erase(runner);
-                                        flag = true;
+                                
+                                // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                                {
+          vector SideA, SideB, SideC, SideH;
+          SideA.CopyVector(&left->second.second->x);
+          SideA.ProjectOntoPlane(&AxisVector);
+//          *out << "SideA: ";
+//          SideA.Output(out);
+//          *out << endl;
+          
+          SideB.CopyVector(&right->second.second->x);
+          SideB.ProjectOntoPlane(&AxisVector);
+//          *out << "SideB: ";
+//          SideB.Output(out);
+//          *out << endl;
+          
+          SideC.CopyVector(&left->second.second->x);
+          SideC.SubtractVector(&right->second.second->x);
+          SideC.ProjectOntoPlane(&AxisVector);
+//          *out << "SideC: ";
+//          SideC.Output(out);
+//          *out << endl;
+
+          SideH.CopyVector(&runner->second.second->x);
+          SideH.ProjectOntoPlane(&AxisVector);
+//          *out << "SideH: ";
+//          SideH.Output(out);
+//          *out << endl;
+          
+          // calculate each length
+          double a = SideA.Norm();
+          //double b = SideB.Norm();
+          //double c = SideC.Norm();
+          double h = SideH.Norm();
+          // calculate the angles
+          double alpha = SideA.Angle(&SideH);
+                                double beta = SideA.Angle(&SideC); 
+                                double gamma = SideB.Angle(&SideH);
+                                double delta = SideC.Angle(&SideH);
+                                double MinDistance = a * sin(beta)/(sin(delta)) * (((alpha < M_PI/2.) || (gamma < M_PI/2.)) ? 1. : -1.);
+//                              *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl; 
+                                *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+                                if ((fabs(h/fabs(h) - MinDistance/fabs(MinDistance)) < MYEPSILON) && (h <  MinDistance)) {
+                                        // throw out point
+                                  *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+                                        BoundaryPoints[axis].erase(runner);
+                                        flag = true;
+                                }
                                 }
                         }
@@ -1472 +1567 @@
                 }
         }
-        int BoundaryPointCount = 0;
-        Walker = start;
-        while (Walker->next != end) {
-                Walker = Walker->next;
-                if (AtomList[Walker->nr] < NDIM) {      // enter all neighbouring points
-                        BoundaryPointCount++;
-                }
-        } 
-        
+
         // 3e. Points no more in the total list, have to be thrown out of each axis lists too!
         for (int axis=0; axis<NDIM; axis++)  {
                 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-                        if (AtomList[runner->second.second->nr] < NDIM)
+                        if (AtomList[runner->second.second->nr] < 1) {
+                          *out << Verbose(1) << "Throwing out " << *runner->second.second << " in axial projection of axis " << axis << "." << endl;
                                 BoundaryPoints[axis].erase(runner);
+                        }
                 }
         }
+        
+        // listing boundary points
+        *out << Verbose(1) << "The following atoms are on the boundary:";
+  int BoundaryPointCount = 0;
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if (AtomList[Walker->nr] > 0) {
+      BoundaryPointCount ++;
+      *out << " " << Walker->Name;
+    }
+  }
+  *out << endl;
 
         *out << Verbose(2) << "I found " << BoundaryPointCount << " points on the convex boundary." << endl;
@@ -1493 +1595 @@
         
         // 4. Create each possible line, number them and put them into a list (3 * AtomCount possible)
-        struct lines {
-                atom *x[2];
-                int nr;
-        } LinesList[3 * AtomCount];
-        
-        // initialise reference storage (needed lateron)
-        int **AtomLines = new int *[AtomCount];
-        Walker = start;
-        while (Walker->next != end) {   // go through all points
-                Walker = Walker->next;
-                if (AtomList[Walker->nr] == NDIM) {     // if this is a boundary point
-                        AtomLines[Walker->nr] = new int[NDIM];
-                        for(int axis=0;axis<NDIM;axis++) {
-                                AtomLines[Walker->nr][2*axis+0] = -1;
-                                AtomLines[Walker->nr][2*axis+1] = -1;
-                        }
-                }
-        }
+        LinesMap AllTriangleLines;
         
         // then create each line
+  *out << Verbose(1) << "Creating lines between adjacent boundary points." << endl;
+  atom *Runner = NULL;
   int LineCount = 0;
-        for(int axis=0;axis<NDIM;axis++) {      // go through every axis
-          for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-            // get left and right neighbours
-            Boundaries::iterator NN[2] = runner;
-            if (NN[0] == BoundaryPoints[axis].begin())
-              NN[0] = BoundaryPoints[axis].end();
-            NN[0]--;
-            NN[1] = runner;
-      NN[1]++;
-      if (NN[1] == BoundaryPoints[axis].end())
-        NN[1] = BoundaryPoints[axis].begin();
-      // check those neighbours for their atomic index
-      for (int k=0;k<2;k++) {
-                        if (runner->second.second->nr > NN[k]->second.second->nr) { // check left neighbour on this projected axis
-                                LinesList[LineCount].x[0] = runner->second.second;
-                                LinesList[LineCount].x[1] = NN[k]->second.second;
-                                if (AtomLines[LinesList[LineCount].x[0]->nr][axis] == -1) {
-                                        AtomLines[LinesList[LineCount].x[0]->nr][axis] = LineCount;
-                                } else if (AtomLines[LinesList[LineCount].x[1]->nr][axis] == -1) {
-                                        AtomLines[LinesList[LineCount].x[1]->nr][axis] = LineCount;
-                                } else { 
-                                        *out << Verbose(2) << "Could not store line number " << LineCount << "." << endl;
-                                }
-                                LineCount++;
-                        }
-      }
-                }
-        }
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if (AtomList[Walker->nr] > 0) { // boundary point!
+      *out << Verbose(1) << "Current Walker is " << *Walker << "." << endl;
+      // make a list of all neighbours
+      DistanceMap Distances;
+      double distance = 0;
+      Runner = start;
+      while (Runner->next != end) {
+        Runner = Runner->next;
+        if ((AtomList[Runner->nr] > 0) && (Runner != Walker)) { // don't mess with ourselves!
+          distance = Walker->x.Distance(&Runner->x);    // note this is squared distance
+          Distances.insert ( DistanceNrPair( distance, Runner) );
+          *out << Verbose(2) << "Inserted distance " << distance << " to " << *Runner << " successfully." << endl;
+        }
+      } // end of distance filling while loop
+      // take 3NNs to draw a line in between
+      int Counter = 0;  // counts up to three inserted lines
+      LinesTestPair LinesTest;
+      for (DistanceMap::iterator spanner = Distances.begin(); spanner != Distances.end(); spanner++) {
+        *out << Verbose(1) << "Current distance to insert is " << spanner->first << " to " << spanner->second+1 << "." << endl;
+        LinesTest = AllTriangleLines.insert ( LinesPair( (Walker->nr*AtomCount)+spanner->second->nr, LineCount) );
+        if (!LinesTest.second) {
+          *out << Verbose(2) << "Insertion of line between " << Walker->nr+1 << " and " << spanner->second->nr+1 << " failed, already present line nr. " << LinesTest.first->second << "." << endl; 
+        } else {
+          LineCount++;
+          Counter++;
+          *out << Verbose(2) << "Insertion of line between " << (LinesTest.first->first/AtomCount)+1 << " and " << (LinesTest.first->first%AtomCount)+1 << " successful, inserting mirrored line and increased counter to " << Counter << "." << endl; 
+          LinesTest = AllTriangleLines.insert ( LinesPair( Walker->nr+(spanner->second->nr*AtomCount), LineCount) );
+          if (!LinesTest.second) {
+            *out << Verbose(2) << "Insertion of mirror line between " << Walker->nr << " and " << spanner->second->nr << " failed, already present line nr. " << LinesTest.first->second << "." << endl;
+            *out << Verbose(0) << "SERIOUS ERROR!!!" << endl;
+          }
+        }
+        if (Counter >= 3) { // stop after three lines (if not run out of possible NNs already)
+          *out << Verbose(2) << "Ok, three lines inserted from current walker " << *Walker << ", stopping." << endl;
+          break;
+        }
+      }
+    } // end of if boundary pointy
+  } // end of loop through all atoms
+
+        // listing created lines
+        *out << Verbose(2) << "The following lines were created:";
+        for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++)
+          if ((liner->first/AtomCount) < (liner->first%AtomCount))
+            *out << " " << (liner->first/AtomCount)+1 << "<->" << (liner->first%AtomCount)+1; 
+        *out << endl;
         *out << Verbose(2) << "I created " << LineCount << " edges." << endl;
         
         // 5. Create every possible triangle from the lines (numbering of each line must be i < j < k)
         struct triangles {
-                atom *x1;
-                atom *x2;
-                atom *x3;
+                atom *x[3];
                 int nr;
         } TriangleList[2 * 3 * AtomCount];      // each line can be part in at most 2 triangles
+
+  *out << Verbose(1) << "Creating triangles out of these lines." << endl;
         int TriangleCount = 0;
-        for (int i=0; i<LineCount; i++) {       // go through every line
+  LinetoAtomPair InsertionPair = LinetoAtomPair(-1, NULL);
+  int endpoints[3];
+  for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++) { // go through every line
                 LinetoAtomMap LinetoAtom;
                 // we have two points, check the lines at either end, whether they share a common endpoint
-                for (int endpoint=0;endpoint<2;endpoint++) {
-                        Walker = LinesList[LineCount].x[endpoint];              
-                        for (int axis=0;axis<NDIM;axis++) {
-                                LinetoAtomTestPair LinetoAtomTest;
-                                LinetoAtomPair InsertionPair = LinetoAtomPair(-1, NULL);
-                                struct lines *TempLine = &LinesList[ AtomLines[Walker->nr][axis] ];
-                                if (AtomLines[Walker->nr][axis] != LineCount) {
-                                        if (TempLine->x[0] != Walker) {
-                                                InsertionPair = LinetoAtomPair(TempLine->x[0]->nr, TempLine->x[0]);
-                                        } else if (TempLine->x[1] != Walker) {
-                                                InsertionPair = LinetoAtomPair(TempLine->x[1]->nr, TempLine->x[0]);
-                                        } else {
-                                                *out << Verbose(2) << "Atom " << *Walker << "is both the end of line " << AtomLines[Walker->nr][axis] << "." << endl;
-                                        }
-                                        if (InsertionPair.first != -1)          // insert if present
-                                                LinetoAtomTest = LinetoAtom.insert( InsertionPair );
-                                        if (!LinetoAtomTest.second) {
-                                                // atom is already present in list, hence we have found a possible triangle
-                                                if ((Walker->nr < LinetoAtomTest.first->first) && (LinetoAtomTest.first->first < InsertionPair.first)) {        // check if numbering is in correct order for uniqueness, if so add
-                                                        TriangleList[TriangleCount].x1 = Walker;
-                                                        TriangleList[TriangleCount].x2 = LinetoAtomTest.first->second;
-                                                        TriangleList[TriangleCount].x3 = InsertionPair.second;
-                                                        TriangleCount++;
-                                                } 
-                                        }       //else {   // check if insertion was successful
-          // successfully inserted, hence nothing found yet
-        //} 
-                                }
-                        }
-                }
+    *out << Verbose(1) << "Examining line between " << (liner->first/AtomCount)+1 << " and " << (liner->first%AtomCount)+1 << "." << endl;
+    int enden[3][2];
+    enden[0][0] = (liner->first/AtomCount);
+    enden[0][1] = (liner->first%AtomCount);
+    if (enden[0][0] < enden[0][1]) {
+                for (int endpoint=0;endpoint<2;endpoint++) {
+                        endpoints[0] = enden[0][endpoint];
+                        *out << Verbose(2) << "Current first endpoint is " << endpoints[0]+1 << "." << endl;
+                        
+                        LinesMap::iterator LineRangeStart = AllTriangleLines.lower_bound(endpoints[0]*AtomCount);
+                        for (LinesMap::iterator coach = LineRangeStart; (coach->first/AtomCount) == endpoints[0]; coach++) { // look at all line's other endpoint
+                    enden[1][0] = (coach->first/AtomCount);
+                    enden[1][1] = (coach->first%AtomCount);
+          //*out << Verbose(3) << "Line is " << coach->first << " with nr. " << coach->second << ": " << enden[1][0+1] << "<->" << enden[1][1]+1 << "." << endl;
+                          endpoints[1] = (enden[1][0] != endpoints[0] ) ? enden[1][0] : enden[1][1];
+                          if ((endpoints[1] > endpoints[0]) && (endpoints[1] != enden[0][(endpoint+1)%2])) {  // don't go back the very same line! 
+                          *out << Verbose(3) << "Current second endpoint is " << endpoints[1]+1 << "." << endl;  
+                          
+              LinesMap::iterator LineRangeStart2 = AllTriangleLines.lower_bound(endpoints[1]*AtomCount);
+              for (LinesMap::iterator coacher = LineRangeStart2; (coacher->first/AtomCount) == endpoints[1]; coacher++) { // look at all line's other endpoint
+                enden[2][0] = (coacher->first/AtomCount);
+                enden[2][1] = (coacher->first%AtomCount);
+              //*out << Verbose(3) << "Line is " << coacher->first << " with nr. " << coacher->second << ": " << enden[2][0]+1 << "<->" << enden[2][1]+1 << "." << endl;
+                endpoints[2] = (enden[2][0] != endpoints[1]) ? enden[2][0] : enden[2][1];
+                if ((endpoints[2] > endpoints[1]) && (endpoints[2] != endpoints[0])) {  // don't go back the very same line! 
+                *out << Verbose(4) << "Current third endpoint is " << endpoints[2]+1 << "." << endl;
+                
+                if (endpoints[2] == enden[0][(endpoint+1)%2]) { // jo, closing triangle
+                  *out << Verbose(3) << "MATCH: Triangle is ";
+                  for (int k=0;k<3;k++) {
+                    *out << endpoints[k]+1 << "\t";
+                    TriangleList[TriangleCount].x[k] = FindAtom(endpoints[k]); 
+                  }
+                  *out << endl;
+                  TriangleList[TriangleCount].nr = TriangleCount;
+                  TriangleCount++;
+                } else {
+                  *out << Verbose(3) << "No match!" << endl;
+                }
+              }
+              }
+                        }
+                        }
+                }
+    }
         }
+  // listing created lines
+  *out << Verbose(2) << "The following triangles were created:";
+  for (int i=0;i<TriangleCount;i++) {
+    *out << " " << TriangleList[i].x[0]->Name << "<->" << TriangleList[i].x[1]->Name << "<->" << TriangleList[i].x[2]->Name; 
+  }
+  *out << endl;
         *out << Verbose(2) << "I created " << TriangleCount << " triangles." << endl;
 
         // 6. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+  *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
         double volume = 0.;
         double PyramidVolume = 0.;
@@ -1595 +1727 @@
         double a,b,c;
         for (int i=0; i<TriangleCount; i++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
-                x.CopyVector(&TriangleList[TriangleCount].x1->x);
-                x.SubtractVector(&TriangleList[TriangleCount].x2->x);
-                y.CopyVector(&TriangleList[TriangleCount].x1->x);
-                y.SubtractVector(&TriangleList[TriangleCount].x3->x);
-                a = sqrt(TriangleList[TriangleCount].x1->x.Distance(&TriangleList[TriangleCount].x2->x));
-                b = sqrt(TriangleList[TriangleCount].x1->x.Distance(&TriangleList[TriangleCount].x3->x));
-                c = sqrt(TriangleList[TriangleCount].x3->x.Distance(&TriangleList[TriangleCount].x2->x));
-                G =  c * (sin(x.Angle(&y))*a); // area of tesselated triangle
-                x.MakeNormalVector(&TriangleList[TriangleCount].x1->x, &TriangleList[TriangleCount].x2->x, &TriangleList[TriangleCount].x3->x);
-                y.CopyVector(&x);
-                x.Scale(CenterOfGravity->Projection(&x));
+                x.CopyVector(&TriangleList[i].x[0]->x);
+                x.SubtractVector(&TriangleList[i].x[1]->x);
+                y.CopyVector(&TriangleList[i].x[0]->x);
+                y.SubtractVector(&TriangleList[i].x[2]->x);
+                a = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[1]->x));
+                b = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[2]->x));
+                c = sqrt(TriangleList[i].x[2]->x.Distance(&TriangleList[i].x[1]->x));
+                G =  sqrt( ( (a*a+b*b+c*c)*(a*a+b*b+c*c) - 2*(a*a*a*a + b*b*b*b + c*c*c*c) )/16.); // area of tesselated triangle
+                x.MakeNormalVector(&TriangleList[i].x[0]->x, &TriangleList[i].x[1]->x, &TriangleList[i].x[2]->x);
+                x.Scale(TriangleList[i].x[1]->x.Projection(&x));
                 h = x.Norm(); // distance of CoG to triangle
                 PyramidVolume = (1./3.) * G * h;                // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
+                *out << Verbose(2) << "Area of triangle is " << G << "A^2, height is " << h << " and the volume is " << PyramidVolume << "." << endl;
                 volume += PyramidVolume; 
         }
         *out << Verbose(1) << "The summed volume is " << volume << " " << (IsAngstroem ? "A^3" : "a_0^3") << "." << endl;
 
-        // free reference lists
+  // 7. translate all points back from CoG
+  *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
+        CenterOfGravity->Scale(-1);
+  Walker=start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    Walker->x.Translate(CenterOfGravity);
+  }
+
+  // free reference lists
         delete[](AtomList);
-  for(int i=0;i<AtomCount;i++)
-    delete[](AtomLines[i]);
-  delete[](AtomLines);
 
         return volume;

src/molecules.hpp

@@ -40 +40 @@
 #define KeySet set<int>
 #define NumberValuePair pair<int, double>
-#define Graph map<KeySet, NumberValuePair, KeyCompare >
-#define GraphPair pair<KeySet, NumberValuePair >
+#define Graph map <KeySet, NumberValuePair, KeyCompare >
+#define GraphPair pair <KeySet, NumberValuePair >
 #define KeySetTestPair pair<KeySet::iterator, bool>
 #define GraphTestPair pair<Graph::iterator, bool>
 
-#define DistanceNrPair pair< double, atom* >
+#define DistanceNrPair pair < double, atom* >
+#define DistanceMap multimap < double, atom* >
+#define DistanceTestPair pair < DistanceMap::iterator, bool>
+
 #define Boundaries map <double, DistanceNrPair >
 #define BoundariesPair pair<double, DistanceNrPair >
-#define BoundariesTestPair pair<Boundaries::iterator, bool>
+#define BoundariesTestPair pair< Boundaries::iterator, bool>
+
+#define LinesMap map <int, int>
+#define LinesPair pair <int, int>
+#define LinesTestPair pair < LinesMap::iterator, bool>
+
 #define LinetoAtomMap map < int, atom * > 
 #define LinetoAtomPair pair < int, atom * > 
@@ -404 +412 @@
     char *configpath;
     char *configname;
+    bool FastParsing;
     
     private: