Changeset a8bcea6 for src

Timestamp:

Dec 4, 2008, 3:15:00 PM (16 years ago)

Author:

Christian Neuen <neuen@…>

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:

f714979

Parents:

f683fe

Message:

several changes, now output is created, quality unknown

Location:

src

Files:

• border.cpp (deleted)
• boundary.cpp (modified) (19 diffs)
• boundary.hpp (modified) (4 diffs)
• builder.cpp (modified) (65 diffs)

src/boundary.cpp ¶

@@ -456 +456 @@
 
 
+/*
+ * This function creates the tecplot file, displaying the tesselation of the hull.
+ * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ */
+void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol)
+{
+  // 8. Store triangles in tecplot file
+  if (tecplot != NULL) {
+    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+    *tecplot << "ZONE T=\"TRIANGLES\", N=" <<  TesselStruct->PointsOnBoundaryCount << ", E=" <<  TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+    int *LookupList = new int[mol->AtomCount];
+    for (int i=0;i<mol->AtomCount;i++)
+      LookupList[i] = -1;
+
+    // print atom coordinates
+    *out << Verbose(2) << "The following triangles were created:";
+    int Counter = 1;
+    atom *Walker = NULL;
+    for (PointMap::iterator target =  TesselStruct->PointsOnBoundary.begin(); target !=  TesselStruct->PointsOnBoundary.end(); target++) {
+      Walker = target->second->node;
+      LookupList[Walker->nr] = Counter++;
+      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+    }
+    *tecplot << endl;
+      // print connectivity
+    for (TriangleMap::iterator runner =  TesselStruct->TrianglesOnBoundary.begin(); runner !=  TesselStruct->TrianglesOnBoundary.end(); runner++) {
+      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
+      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+    }
+    delete[](LookupList);
+    *out << endl;
+  }
+}
+
 /** Determines the volume of a cluster.
  * Determines first the convex envelope, then tesselates it and calculates its volume.
@@ -478 +514 @@
   double a,b,c;
 
-  Find_non_convex_border(*TesselStruct, mol);
+  Find_non_convex_border(out, tecplot, *TesselStruct, mol);
   /*
   // 1. calculate center of gravity
@@ -560 +596 @@
 
 
-
-  // 8. Store triangles in tecplot file
-  if (tecplot != NULL) {
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES\", N=" <<  TesselStruct->PointsOnBoundaryCount << ", E=" <<  TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int *LookupList = new int[mol->AtomCount];
-    for (int i=0;i<mol->AtomCount;i++)
-      LookupList[i] = -1;
-
-    // print atom coordinates
-    *out << Verbose(2) << "The following triangles were created:";
-    int Counter = 1;
-    atom *Walker = NULL;
-    for (PointMap::iterator target =  TesselStruct->PointsOnBoundary.begin(); target !=  TesselStruct->PointsOnBoundary.end(); target++) {
-      Walker = target->second->node;
-      LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
-    }
-    *tecplot << endl;
-      // print connectivity
-    for (TriangleMap::iterator runner =  TesselStruct->TrianglesOnBoundary.begin(); runner !=  TesselStruct->TrianglesOnBoundary.end(); runner++) {
-      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
-    }
-    delete[](LookupList);
-    *out << endl;
-  }
+  write_tecplot_file(out, tecplot, TesselStruct, mol);
+
 
   // free reference lists
@@ -699 +709 @@
 Tesselation::~Tesselation()
 {
+        cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
   for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
     delete(runner->second);
@@ -1028 +1039 @@
  *  with all neighbours of the candidate.
  */
-void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule* mol)
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
 {
   /* d2 is normal vector on the triangle
@@ -1045 +1056 @@
   AngleCheck.ProjectOntoPlane(Chord);
 
-  if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  if (problem)
   {
-
-    if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
+          cout << "Atom number" << Candidate->nr << endl;
+          Candidate->x.Output((ofstream *)&cout);
+          cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
+  }
+
+  if (a != Candidate and b != Candidate)
+  {
+
+          if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+          {
+                  if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[0]) //This will give absolute preference to those in "right-hand" quadrants
+                  {
+                          Opt_Candidate = Candidate;
+                          Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+                          Storage[1]=AngleCheck.Angle(d2);
+                  }
+                  else
+                  {
+                          if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]< AngleCheck.Angle(d2)) or \
+                                          (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]> AngleCheck.Angle(d2)))
+                                  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+                          {
+                                  Opt_Candidate = Candidate;
+                                  Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+                                  Storage[1]=AngleCheck.Angle(d2);
+                          }
+                          else
+                          {
+                                  if (problem)
+                                          cout << "Loses to better candidate" << endl;
+                          }
+                  }
+          }
+          else
+          {
+                  if (problem)
+                          cout << "erfuellt Dreiecksbedingung fuer sehne nicht" <<endl;
+          }
+  }
+  else
+  {
+          if (problem)
+                          cout << "identisch mit Ursprungslinie" << endl;
+  }
+
+  if (n<5) // Five is the recursion level threshold.
+  {
+          for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
       {
-        Storage[0]=(double)Candidate->nr;
-        Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-        Storage[2]=AngleCheck.Angle(d2);
+                  Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+
+                  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
       }
-    else
-      {
-        if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
-            (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
-          //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-          {
-            Storage[0]=(double)Candidate->nr;
-            Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-            Storage[2]=AngleCheck.Angle(d2);
-          }
-      }
-  }
-
-  if (n<5) // Five is the recursion level threshold.
-    {
-      for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
-        {
-                Walker= mol->start; // go through all atoms
-
-          while (Walker->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr))
-            { // until atom found which belongs to bond
-              Walker = Walker->next;
-            }
-
-
-          Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Storage, RADIUS, mol);  //call function again
-        }
-    }
+  }
 };
 
@@ -1087 +1117 @@
  * this function fins a triangle to a line, adjacent to an existing one.
  */
-void Tesselation::Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
-{
-        printf("Looking for next suitable triangle \n");
+void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS)
+{
+        cout << Verbose(1) << "Looking for next suitable triangle \n";
   Vector direction1;
   Vector helper;
@@ -1095 +1125 @@
   atom* Walker;
 
-  double *Storage;
-  Storage = new double[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
-
-
+  double Storage[2];
+  Storage[0]=-2.;    // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[1]=9999999.;  // This is also lower then any value produced by an eligible atom, which are all positive
+  atom* Opt_Candidate = NULL;
+
+
+  cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
   helper.CopyVector(&(Line.endpoints[0]->node->x));
   for (int i =0; i<3; i++)
@@ -1115 +1145 @@
   direction1.CopyVector(&Line.endpoints[0]->node->x);
   direction1.SubtractVector(&Line.endpoints[1]->node->x);
-  direction1.VectorProduct(T.NormalVector);
+  direction1.VectorProduct(&(T.NormalVector));
 
   if (direction1.ScalarProduct(&helper)>0)
@@ -1125 +1155 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-printf("Looking for third point candidates for triangle \n");
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, T.NormalVector, Storage, RADIUS, mol);
-
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
+
+  if (Opt_Candidate==NULL)
+  {
+          cout << Verbose(1) << "No new Atom found, triangle construction will crash" << endl;
+          write_tecplot_file(out, tecplot, this, mol);
+          tecplot->flush();
+          Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol, 1);
+
+  }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
-  // Next Triangle is Line, atom with number in Storage[0]
-
-  Walker= mol->start;
-  while (Walker->nr != (int)Storage[0])
-    {
-      Walker = Walker->next;
-    }
-
-  AddPoint(Walker);
-
-  BPS[0] = new class BoundaryPointSet(Walker);
+
+  cout << Verbose(1) << "Adding exactly one Walker for reasons completely unknown to me ... " << endl;
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  AddPoint(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(Walker);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
@@ -1151 +1189 @@
   TrianglesOnBoundaryCount++;
 
+  cout << Verbose(1) << "Checking for already present lines ... " << endl;
   for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
     {
@@ -1160 +1199 @@
         }
     }
-  BTS->GetNormalVector(*BTS->NormalVector);
-
-  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
-      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
+  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
+
+  BTS->GetNormalVector(BTS->NormalVector);
+
+  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector))<0 && Storage[0]>0) ||
+                  ((BTS->NormalVector.ScalarProduct(&(T.NormalVector))>0 && Storage[0]<0)) )
     {
-      BTS->NormalVector->Scale(-1);
-    }
-
-};
-
-
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, double Storage[3], molecule* mol)
-{
-        printf("Looking for second point of starting triangle \n");
-  int i;
-  Vector *AngleCheck;
-  atom* Walker;
-
-  AngleCheck->CopyVector(&(Candidate->x));
-  AngleCheck->SubtractVector(&(a->x));
-  if (AngleCheck->Angle(&Oben) < Storage[1])
+      BTS->NormalVector.Scale(-1);
+    };
+
+};
+
+
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, atom*& Opt_Candidate, double Storage[2], molecule* mol, double RADIUS)
+{
+        cout << Verbose(1) << "Looking for second point of starting triangle, recursive level "<< n <<endl;;
+        int i;
+        Vector AngleCheck;
+        atom* Walker;
+
+        AngleCheck.CopyVector(&(Candidate->x));
+        AngleCheck.SubtractVector(&(a->x));
+        if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
     {
-          //printf("Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-      Storage[0]=(double)(Candidate->nr);
-      Storage[1]=AngleCheck->Angle(&Oben);
-      //printf("Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+          //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+      Opt_Candidate=Candidate;
+      Storage[0]=AngleCheck.Angle(&Oben);
+      //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
     };
-  printf("%d \n", n);
-  if (n<5)
+        if (n<5)
     {
       for (i = 0; i< mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
-        {
-          Walker = mol->start;
-          while (Candidate->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr))
-            {
-              Walker = Walker->next;
-            };
-
-          Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Storage, mol);
-            };
+      {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+          Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+      };
     };
 
@@ -1208 +1243 @@
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-        printf("Looking for starting triangle \n");
+        cout << Verbose(1) << "Looking for starting triangle \n";
   int i=0;
   atom* Walker;
@@ -1227 +1262 @@
       max_coordinate[i] =-1;
     }
-printf("Molecule mol is there and has %d Atoms \n", mol->AtomCount);
+cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
   Walker = mol->start;
   while (Walker->next != mol->end)
@@ -1242 +1277 @@
     }
 
-printf("Found starting atom \n");
+cout << Verbose(1) << "Found starting atom \n";
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
   const int k=2;
@@ -1253 +1288 @@
       Walker = Walker->next;
     }
-printf("%d \n", Walker->nr);
-  double Storage[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[2]=999999.;  // This will be an angle looking for the third point.
-printf("%d \n", mol->NumberOfBondsPerAtom[Walker->nr]);
+cout << Verbose(1) << "Number of start atom: " << Walker->nr << endl;
+  double Storage[2];
+  atom* Opt_Candidate = NULL;
+  Storage[0]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[1]=999999.;  // This will be an angle looking for the third point.
+cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[Walker->nr] << endl;
 
   for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
     {
-      Walker2 = mol->start;
-      Walker2 = Walker2->next;
-      while (Walker2->nr != (mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr) )
-        {
-          Walker2 = Walker2->next;
-        }
-
-      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
-    }
-
-  Walker2 = mol->start;
-
-  while (Walker2->nr != int(Storage[0]))
-    {
-      Walker2 = Walker2->next;
-    }
-
+          Walker2 = mol->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
+    }
+
+  Walker2 = Opt_Candidate;
+  Opt_Candidate=NULL;
   helper.CopyVector(&(Walker->x));
   helper.SubtractVector(&(Walker2->x));
   Oben.ProjectOntoPlane(&helper);
   helper.VectorProduct(&Oben);
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[2]= -10.;  // This will be an angle looking for the third point.
+  Storage[0]=-2.;       // This will indicate the quadrant.
+  Storage[1]= 9999999.; // This will be an angle looking for the third point.
 
   Chord.CopyVector(&(Walker->x)); // bring into calling function
   Chord.SubtractVector(&(Walker2->x));
 
-printf("Looking for third point candidates \n");
-       Find_next_suitable_point(Walker, Walker2, (mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][0]->rightatom : mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom), 0, &Chord, &helper, &Oben, Storage, RADIUS, mol);
-  Walker3 = mol->start;
-  while (Walker3->nr != int(Storage[0]))
-    {
-      Walker3 = Walker3->next;
-    }
-
-  //Starting Triangle is Walker, Walker2, index Storage[0]
+  cout << Verbose(1) << "Looking for third point candidates \n";
+  Find_next_suitable_point(Walker, Walker2, mol->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
+
+
+  //Starting Triangle is Walker, Walker2, Opt_Candidate
 
   AddPoint(Walker);
   AddPoint(Walker2);
-  AddPoint(Walker3);
+  AddPoint(Opt_Candidate);
+
+  cout << Verbose(1) << "The found starting triangle consists of " << *Walker << ", " << *Walker2 << " and " << *Opt_Candidate << "." << endl;
 
   BPS[0] = new class BoundaryPointSet(Walker);
@@ -1308 +1329 @@
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
-  BPS[1] = new class BoundaryPointSet(Walker3);
+  BPS[1] = new class BoundaryPointSet(Opt_Candidate);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
@@ -1314 +1335 @@
   BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  Tesselation::BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
   TrianglesOnBoundaryCount++;
@@ -1324 +1345 @@
     };
 
-       BTS->GetNormalVector(*BTS->NormalVector);
-
-       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
+       BTS->GetNormalVector(BTS->NormalVector);
+
+       if( BTS->NormalVector.ScalarProduct(&Oben)<0)
          {
-           BTS->NormalVector->Scale(-1);
+           BTS->NormalVector.Scale(-1);
          }
 };
 
 
-void Find_non_convex_border(Tesselation Tess, molecule* mol)
-{
-        printf("Entering finding of non convex hull. \n");
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol)
+{
+        cout << Verbose(1) << "Entering finding of non convex hull. " << endl;
+        cout << flush;
   const double RADIUS =6;
+  LineMap::iterator baseline;
   Tess.Find_starting_triangle(mol, RADIUS);
 
-  for (LineMap::iterator baseline = Tess.LinesOnBoundary.begin(); baseline != Tess.LinesOnBoundary.end(); baseline++)
+  baseline = Tess.LinesOnBoundary.begin();
+  while (baseline != Tess.LinesOnBoundary.end()) {
     if (baseline->second->TrianglesCount == 1)
       {
-                Tess.Find_next_suitable_triangle(mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+        cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+                Tess.Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+        cout << Verbose(1) << "End of Tesselation ... " << endl;
       }
     else
       {
-        printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+        cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
       }
-};
+  baseline++;
+  }
+
+};

src/boundary.hpp ¶

@@ -36 +36 @@
     BoundaryPointSet(atom *Walker);
     ~BoundaryPointSet();
-    
+
     void AddLine(class BoundaryLineSet *line);
-    
+
     LineMap lines;
     int LinesCount;
@@ -64 +64 @@
     BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
     ~BoundaryTriangleSet();
-    
+
     void GetNormalVector(Vector &NormalVector);
-    
+
     class BoundaryPointSet *endpoints[3];
     class BoundaryLineSet *lines[3];
-    Vector *NormalVector;
+    Vector NormalVector;
     int Nr;
 };
@@ -75 +75 @@
 class Tesselation {
   public:
-    
+
     Tesselation();
     ~Tesselation();
-    
+
     void TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol);
     void GuessStartingTriangle(ofstream *out);
     void AddPoint(atom * Walker);
     void Find_starting_triangle(molecule* mol, const double RADIUS);
-    void Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS);
-    
+    void Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS);
+
     PointMap PointsOnBoundary;
     LineMap LinesOnBoundary;
@@ -105 +105 @@
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule *mol);
-void Find_non_convex_border(Tesselation Tess, molecule* mol);
+void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule *mol, bool problem);
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol);
 
 

src/builder.cpp ¶

@@ -1 +1 @@
 /** \file builder.cpp
- * 
+ *
  * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
  * The output is the complete configuration file for PCP for direct use.
@@ -6 +6 @@
  * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
  * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
- *    
+ *
  */
 
 /*! \mainpage Molecuilder - a molecular set builder
- * 
+ *
  * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
- * 
+ *
  * \section about About the Program
- * 
+ *
  *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
  *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
  *  already constructed atoms.
- * 
+ *
  *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
  *  molecular dynamics implementation.
- * 
+ *
  * \section install Installation
- * 
+ *
  *  Installation should without problems succeed as follows:
  *  -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
  *  -# make
  *  -# make install
- * 
+ *
  *  Further useful commands are
  *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
- *  -# make doxygen-doc: Creates these html pages out of the documented source 
- * 
+ *  -# make doxygen-doc: Creates these html pages out of the documented source
+ *
  * \section run Running
- * 
+ *
  *  The program can be executed by running: ./molecuilder
- * 
+ *
  *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
  *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
- *  later re-execution. 
- * 
+ *  later re-execution.
+ *
  * \section ref References
- *  
+ *
  *  For the special configuration file format, see the documentation of pcp.
- * 
+ *
  */
 
@@ -80 +80 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
       case 'a': // absolute coordinates of atom
@@ -89 +89 @@
         mol->AddAtom(first);  // add to molecule
         break;
-        
+
       case 'b': // relative coordinates of atom wrt to reference point
         first = new atom;
@@ -105 +105 @@
         mol->AddAtom(first);  // add to molecule
         break;
-        
+
       case 'c': // relative coordinates of atom wrt to already placed atom
         first = new atom;
@@ -111 +111 @@
         do {
           if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
-          second = mol->AskAtom("Enter atom number: ");                
+          second = mol->AskAtom("Enter atom number: ");
           cout << Verbose(0) << "Enter relative coordinates." << endl;
           first->x.AskPosition(mol->cell_size, false);
@@ -121 +121 @@
         mol->AddAtom(first);  // add to molecule
         break;
-      
+
       case 'd': // two atoms, two angles and a distance
         first = new atom;
@@ -152 +152 @@
           x.Copyvector(&fourth->x);
           x.SubtractVector(&third->x);
-          
+
           if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
             cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
@@ -167 +167 @@
           cout << "x: ",
           x.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           z.MakeNormalVector(&second->x,&third->x,&fourth->x);
           cout << "z: ",
           z.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           y.MakeNormalVector(&x,&z);
           cout << "y: ",
           y.Output((ofstream *)&cout);
-          cout << endl;          
-         
+          cout << endl;
+
           // rotate vector around first angle
           first->x.CopyVector(&x);
@@ -182 +182 @@
           cout << "Rotated vector: ",
           first->x.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           // remove the projection onto the rotation plane of the second angle
           n.CopyVector(&y);
@@ -188 +188 @@
           cout << "N1: ",
           n.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           first->x.SubtractVector(&n);
           cout << "Subtracted vector: ",
           first->x.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           n.CopyVector(&z);
           n.Scale(first->x.Projection(&z));
           cout << "N2: ",
           n.Output((ofstream *)&cout);
-          cout << endl;          
+          cout << endl;
           first->x.SubtractVector(&n);
           cout << "2nd subtracted vector: ",
           first->x.Output((ofstream *)&cout);
-          cout << endl;          
-          
+          cout << endl;
+
           // rotate another vector around second angle
           n.CopyVector(&y);
@@ -208 +208 @@
           cout << "2nd Rotated vector: ",
           n.Output((ofstream *)&cout);
-          cout << endl;          
-          
+          cout << endl;
+
           // add the two linear independent vectors
           first->x.AddVector(&n);
-          first->x.Normalize();         
+          first->x.Normalize();
           first->x.Scale(a);
           first->x.AddVector(&second->x);
-          
+
           cout << Verbose(0) << "resulting coordinates: ";
           first->x.Output((ofstream *)&cout);
@@ -241 +241 @@
         } while ((j != -1) && (i<128));
         if (i >= 2) {
-          first->x.LSQdistance(atoms, i);             
+          first->x.LSQdistance(atoms, i);
 
           first->x.Output((ofstream *)&cout);
@@ -261 +261 @@
   Vector x, y;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
   cout << Verbose(0) << " a - on origin" << endl;
@@ -271 +271 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -327 +327 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -346 +346 @@
       second = mol->AskAtom("Enter second atom: ");
 
-      n.CopyVector((const Vector *)&first->x);  
-      n.SubtractVector((const Vector *)&second->x);  
+      n.CopyVector((const Vector *)&first->x);
+      n.SubtractVector((const Vector *)&second->x);
       n.Normalize();
-      break;        
+      break;
     case 'd':
       char shorthand[4];
@@ -363 +363 @@
         x.x[i] = gsl_vector_get(param.x,i);
         n.x[i] = gsl_vector_get(param.x,i+NDIM);
-      } 
+      }
       gsl_vector_free(param.x);
       cout << Verbose(0) << "Offset vector: ";
@@ -369 +369 @@
       cout << Verbose(0) << endl;
       n.Normalize();
-      break;       
+      break;
   };
   cout << Verbose(0) << "Alignment vector: ";
@@ -385 +385 @@
   Vector n;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
   cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
@@ -394 +394 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -413 +413 @@
       second = mol->AskAtom("Enter second atom: ");
 
-      n.CopyVector((const Vector *)&first->x);  
-      n.SubtractVector((const Vector *)&second->x);  
+      n.CopyVector((const Vector *)&first->x);
+      n.SubtractVector((const Vector *)&second->x);
       n.Normalize();
-      break;          
+      break;
   };
   cout << Verbose(0) << "Normal vector: ";
@@ -433 +433 @@
   double tmp1, tmp2;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
   cout << Verbose(0) << " a - state atom for removal by number" << endl;
@@ -442 +442 @@
   cout << Verbose(0) << "INPUT: ";
   cin >> choice;
-  
+
   switch (choice) {
     default:
@@ -475 +475 @@
           mol->RemoveAtom(first);
       }
-      break;          
+      break;
   };
   //mol->Output((ofstream *)&cout);
@@ -492 +492 @@
   int Z;
   char choice;  // menu choice char
-  
+
   cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
   cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
@@ -514 +514 @@
       for (int i=MAX_ELEMENTS;i--;)
         min[i] = 0.;
-        
-      second = mol->start;   
+
+      second = mol->start;
       while ((second->next != mol->end)) {
         second = second->next; // advance
@@ -526 +526 @@
         }
         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
-        //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;         
+        //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
       }
       for (int i=MAX_ELEMENTS;i--;)
         if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
       break;
-      
+
     case 'b':
       first = mol->AskAtom("Enter first atom: ");
@@ -556 +556 @@
       y.SubtractVector((const Vector *)&second->x);
       cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
-      cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;         
+      cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
       break;
     case 'd':
@@ -600 +600 @@
   int Order1;
   clock_t start, end;
-  
+
   cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
   cout << Verbose(0) << "What's the desired bond order: ";
@@ -609 +609 @@
     end = clock();
     cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
-  } else 
+  } else
     cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
 };
@@ -623 +623 @@
   atom *Walker = mol->start;
   int i, comp, counter=0;
-  
+
   // generate some KeySets
   cout << "Generating KeySets." << endl;
@@ -637 +637 @@
   cout << "Testing insertion of already present item in KeySets." << endl;
   KeySetTestPair test;
-  test = TestSets[mol->AtomCount-1].insert(Walker->nr); 
+  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
   if (test.second) {
     cout << Verbose(1) << "Insertion worked?!" << endl;
@@ -646 +646 @@
   TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
 
-  // constructing Graph structure  
+  // constructing Graph structure
   cout << "Generating Subgraph class." << endl;
   Graph Subgraphs;
@@ -657 +657 @@
   cout << "Testing insertion of already present item in Subgraph." << endl;
   GraphTestPair test2;
-  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 
+  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
   if (test2.second) {
     cout << Verbose(1) << "Insertion worked?!" << endl;
@@ -663 +663 @@
     cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
   }
-  
+
   // show graphs
   cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
@@ -674 +674 @@
       if ((*key) > comp)
         cout << (*key) << " ";
-      else 
+      else
         cout << (*key) << "! ";
       comp = (*key);
@@ -716 +716 @@
   else
     cout << "failed." << endl;
-  
+
   // and save to xyz file
   if (ConfigFileName != NULL) {
@@ -727 +727 @@
     strcat(filename, ".xyz");
     output.open(filename, ios::trunc);
-  } 
+  }
   cout << Verbose(0) << "Saving of XYZ file ";
   if (mol->MDSteps <= 1) {
@@ -742 +742 @@
   output.close();
   output.clear();
-  
+
   // and save as MPQC configuration
   if (ConfigFileName != NULL)
@@ -753 +753 @@
   else
     cout << "failed." << endl;
-  
+
   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
     cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
@@ -785 +785 @@
   int argptr;
   PathToDatabases = LocalPath;
-  
+
   if (argc > 1) { // config file specified as option
     // 1. : Parse options that just set variables or print help
@@ -798 +798 @@
           case '?':
             cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
-            cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl; 
+            cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
             cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
             cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
+            cout << "\t-A <source>\tCreate adjacenzy list from bonds parsed from 'dbond'-style file." <<endl;
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
@@ -812 +813 @@
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
             cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
-            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
+            cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
-            cout << "\t-T <file> Store temperatures from the config file in <file>." << endl; 
+            cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
-            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 
+            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
             cout << "\t-v/-V\t\tGives version information." << endl;
             cout << "Note: config files must not begin with '-' !" << endl;
@@ -854 +855 @@
         argptr++;
     } while (argptr < argc);
-    
+
     // 2. Parse the element database
     if (periode->LoadPeriodentafel(PathToDatabases)) {
@@ -863 +864 @@
       return 1;
     }
-    
+
     // 3. Find config file name and parse if possible
     if (argv[1][0] != '-') {
@@ -902 +903 @@
     } else
       config_present = absent;
-    
+
     // 4. parse again through options, now for those depending on elements db and config presence
     argptr = 1;
@@ -946 +947 @@
                     config_present = present;
                 } else
-                  cerr << Verbose(1) << "Could not find the specified element." << endl; 
+                  cerr << Verbose(1) << "Could not find the specified element." << endl;
                 argptr+=4;
               }
@@ -956 +957 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
-            case 'D': 
+            case 'D':
               ExitFlag = 1;
               {
@@ -1002 +1003 @@
               }
               break;
+            case 'A':
+                ExitFlag = 1;
+                if ((argptr >= argc) || (argv[argptr][0] == '-')){
+                        ExitFlag =255;
+                        cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
+                }
+                else{
+                        cout << "Parsing bonds from " << argv[argptr] << "." << endl;
+                    ifstream *input = new ifstream(argv[argptr]);
+                        mol->CreateAdjacencyList2((ofstream *)&cout, input);
+                        input->close();
+                }
+                break;
             case 'T':
               ExitFlag = 1;
@@ -1168 +1182 @@
             case 'o':
               ExitFlag = 1;
-              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+              if ((argptr >= argc) || (argv[argptr][0] == '-')){
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
               } else {
                 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+                ofstream *output = new ofstream(argv[argptr], ios::trunc);
+                //$$$
                 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
-                ofstream *output = new ofstream(argv[argptr], ios::trunc);
                 VolumeOfConvexEnvelope((ofstream *)&cout, output, &configuration, NULL, mol);
                 output->close();
@@ -1269 +1284 @@
                       mol->Translate(&x);
                     }
-                    mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+                    mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
                   }
                 }
@@ -1332 +1347 @@
   if (j == 1) return 0; // just for -v and -h options
   if (j) return j;  // something went wrong
-  
+
   // General stuff
   if (mol->cell_size[0] == 0.) {
@@ -1346 +1361 @@
   // now the main construction loop
   cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
-  do {    
+  do {
     cout << Verbose(0) << endl << endl;
     cout << Verbose(0) << "============Element list=======================" << endl;
@@ -1365 +1380 @@
     cout << Verbose(0) << "-----------------------------------------------" << endl;
     cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
-    cout << Verbose(0) << "i - realign molecule" << endl; 
-    cout << Verbose(0) << "m - mirror all molecules" << endl; 
+    cout << Verbose(0) << "i - realign molecule" << endl;
+    cout << Verbose(0) << "m - mirror all molecules" << endl;
     cout << Verbose(0) << "t - translate molecule by vector" << endl;
     cout << Verbose(0) << "c - scale by unit transformation" << endl;
@@ -1377 +1392 @@
     cout << Verbose(0) << "Input: ";
     cin >> choice;
-    
+
     switch (choice) {
       default:
       case 'a': // add atom
         AddAtoms(periode, mol);
-        choice = 'a'; 
-        break;
-      
+        choice = 'a';
+        break;
+
       case 'b': // scale a bond
         cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
@@ -1400 +1415 @@
         }
         //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
-        //second->Output(second->type->No, 1, (ofstream *)&cout);       
-        break;
-
-      case 'c': // unit scaling of the metric 
+        //second->Output(second->type->No, 1, (ofstream *)&cout);
+        break;
+
+      case 'c': // unit scaling of the metric
        cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
        cout << Verbose(0) << "Enter three factors: ";
@@ -1414 +1429 @@
        delete[](factor);
        break;
-       
+
       case 'd': // duplicate the periodic cell along a given axis, given times
         cout << Verbose(0) << "State the axis [(+-)123]: ";
@@ -1420 +1435 @@
         cout << Verbose(0) << "State the factor: ";
         cin >> faktor;
-        
+
         mol->CountAtoms((ofstream *)&cout);  // recount atoms
         if (mol->AtomCount != 0) {  // if there is more than none
@@ -1461 +1476 @@
             mol->Translate(&x);
           }
-          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
         }
         break;
-      
+
       case 'e': // edit each field of the configuration
        configuration.Edit(mol);
        break;
-  
+
       case 'f':
         FragmentAtoms(mol, &configuration);
         break;
-        
+
       case 'g': // center the atoms
         CenterAtoms(mol);
         break;
-        
-      case 'i': // align all atoms 
+
+      case 'i': // align all atoms
         AlignAtoms(periode, mol);
         break;
@@ -1488 +1503 @@
         MirrorAtoms(mol);
         break;
-       
+
       case 'o': // create the connection matrix
         {
@@ -1509 +1524 @@
         }
         break;
-        
+
       case 'p': // parse and XYZ file
         cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
@@ -1518 +1533 @@
         break;
 
-      case 'q': // quit 
-        break;
-        
+      case 'q': // quit
+        break;
+
       case 'r': // remove atom
-        RemoveAtoms(mol);        
-        break;
-        
+        RemoveAtoms(mol);
+        break;
+
       case 's': // save to config file
         SaveConfig(ConfigFileName, &configuration, periode, mol);
@@ -1530 +1545 @@
 
       case 't': // translate all atoms
-       cout << Verbose(0) << "Enter translation vector." << endl;       
+       cout << Verbose(0) << "Enter translation vector." << endl;
        x.AskPosition(mol->cell_size,0);
        mol->Translate((const Vector *)&x);
        break;
-  
+
       case 'T':
         testroutine(mol);
         break;
-      
+
       case 'u': // change an atom's element
         first = NULL;
@@ -1545 +1560 @@
           cin >> Z;
         } while ((first = mol->FindAtom(Z)) == NULL);
-        cout << Verbose(0) << "New element by atomic number Z: ";        
+        cout << Verbose(0) << "New element by atomic number Z: ";
         cin >> Z;
         first->type = periode->FindElement(Z);
-        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;    
+        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
         break;
     };
   } while (choice != 'q');
-  
+
   // save element data base
   if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName