Changeset d52ea1b

Timestamp:

Jun 7, 2008, 1:21:53 PM (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:

17ce5c

Parents:

65684f

Message:

Working version of PAS transformation (tested on C-S-H cluster) and code scaffold of VolumeOfConvexEnvelope that's untested so far

PAS simply calculates inertia tensor, diagonalizes it via GSL and transforms molecule into eigenvector system such that z axis is eigenvector of biggest eigenvalue.

Location:

src

Files:

• builder.cpp (modified) (7 diffs)
• molecules.cpp (modified) (24 diffs)
• molecules.hpp (modified) (5 diffs)

src/builder.cpp

@@ -484 +484 @@
  * \param *mol the molecule with all the atoms
  */
-static void MeasureAtoms(periodentafel *periode, molecule *mol)
+static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
 {
   atom *first, *second, *third;
@@ -496 +496 @@
   cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
   cout << Verbose(0) << " c - calculate bond angle" << endl;
+  cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
+  cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
   cout << Verbose(0) << "all else - go back" << endl;
   cout << Verbose(0) << "===============================================" << endl;
@@ -553 +555 @@
       cout << Verbose(0) << (acos(x.ScalarProduct((const vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;         
       break;
+    case 'd':
+        cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+        cout << Verbose(0) << "Shall we rotate? [0/1]: ";
+        cin >> Z;
+        if ((Z >=0) && (Z <=1))
+          mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
+        else
+          mol->PrincipalAxisSystem((ofstream *)&cout, false);
+        break;
+    case 'e':
+        cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+        mol->VolumeOfConvexEnvelope((ofstream *)&cout, configuration->GetIsAngstroem());
+        break;
   }
 };
@@ -745 +760 @@
             cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
+            cout << "\t-m\tAlign in PAS with greatest EV along z axis." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
@@ -945 +961 @@
               }
               break;
+            case 'm':
+              ExitFlag = 1;
+              cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+              mol->PrincipalAxisSystem((ofstream *)&cout, true);
+              break;
             default:   // no match? Step on
               argptr++;
@@ -1151 +1172 @@
 
       case 'l': // measure distances or angles
-        MeasureAtoms(periode, mol);
+        MeasureAtoms(periode, mol, &configuration);
         break;
 
@@ -1164 +1185 @@
         mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
         //mol->CreateListOfBondsPerAtom((ofstream *)&cout);
-        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
+        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, MinimumRingSize);
         while (Subgraphs->next != NULL) {
           Subgraphs = Subgraphs->next;

src/molecules.cpp

@@ -179 +179 @@
   bond *FirstBond = NULL, *SecondBond = NULL; // Other bonds in double bond case to determine "other" plane
   atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
-  int i;  // loop variable
   double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
   vector OrthoVector1, OrthoVector2;  // temporary vectors in coordination construction
@@ -253 +252 @@
     case 2:
       // determine two other bonds (warning if there are more than two other) plus valence sanity check
-      for (i=0;i<NumBond;i++) {
+      for (int i=0;i<NumBond;i++) {
         if (BondList[i] != TopBond) {
           if (FirstBond == NULL) {
@@ -312 +311 @@
       FirstOtherAtom->x.Zero();
       SecondOtherAtom->x.Zero();
-      for(i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
+      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
         FirstOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (sin(bondangle));
         SecondOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (-sin(bondangle));
@@ -319 +318 @@
       SecondOtherAtom->x.Scale(&BondRescale);
       //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
-      for(i=NDIM;i--;) { // and make relative to origin atom
+      for(int i=NDIM;i--;) { // and make relative to origin atom
         FirstOtherAtom->x.x[i] += TopOrigin->x.x[i]; 
         SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
@@ -439 +438 @@
   int NumberOfAtoms = 0; // atom number in xyz read
   int i, j; // loop variables
-  atom *first = NULL;  // pointer to added atom
+  atom *Walker = NULL;  // pointer to added atom
   char shorthand[3];  // shorthand for atom name
   ifstream xyzfile;   // xyz file
@@ -457 +456 @@
 
   for(i=0;i<NumberOfAtoms;i++){
-    first = new atom;
+    Walker = new atom;
     getline(xyzfile,line,'\n');
     istringstream *item = new istringstream(line);
@@ -466 +465 @@
     *item >> x[1];
     *item >> x[2];
-    first->type = elemente->FindElement(shorthand);
-    if (first->type == NULL) {
+    Walker->type = elemente->FindElement(shorthand);
+    if (Walker->type == NULL) {
       cerr << "Could not parse the element at line: '" << line << "', setting to H.";
-      first->type = elemente->FindElement(1);
+      Walker->type = elemente->FindElement(1);
     }
     for(j=NDIM;j--;)
-      first->x.x[j] = x[j];
-     AddAtom(first);  // add to molecule
+      Walker->x.x[j] = x[j];
+     AddAtom(Walker);  // add to molecule
     delete(item);
   }
@@ -710 +709 @@
 }; 
 
-/** Centers the center of gravity of the atoms at (0,0,0).
+/** Returns vector pointing to center of gravity.
  * \param *out output stream for debugging
- * \param *center return vector for translation vector
- */
-void molecule::CenterGravity(ofstream *out, vector *center)
-{
+ * \return pointer to center of gravity vector
+ */
+vector * molecule::DetermineCenterOfGravity(ofstream *out)
+{
+        atom *ptr = start->next;  // start at first in list
+        vector *a = new vector();
+        vector tmp;
   double Num = 0;
-  atom *ptr = start->next;  // start at first in list
-  vector tmp;
-  
-  for(int i=NDIM;i--;) // zero center vector
-    center->x[i] = 0.;
-    
+        
+        a->Zero();
+
   if (ptr != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
@@ -729 +728 @@
       tmp.CopyVector(&ptr->x);
       tmp.Scale(ptr->type->mass);  // scale by mass
-      center->AddVector(&tmp);       
-    }
-    center->Scale(-1./Num); // divide through total mass (and sign for direction)
+      a->AddVector(&tmp);       
+    }
+    a->Scale(-1./Num); // divide through total mass (and sign for direction)
+  }
+  *out << Verbose(1) << "Resulting center of gravity: ";
+  a->Output(out);
+  *out << endl;
+  return a;
+};
+
+/** Centers the center of gravity of the atoms at (0,0,0).
+ * \param *out output stream for debugging
+ * \param *center return vector for translation vector
+ */
+void molecule::CenterGravity(ofstream *out, vector *center)
+{
+  if (center == NULL) {
+    DetermineCenter(*center);
+    Translate(center);
+    delete(center);
+  } else {
     Translate(center);
   }
@@ -775 +792 @@
 };
 
-/** Determines center of gravity (yet not considering atom masses).
- * \param CenterOfGravity reference to return vector
- */
-void molecule::DetermineCenterOfGravity(vector &CenterOfGravity)
+/** Determines center of molecule (yet not considering atom masses).
+ * \param Center reference to return vector
+ */
+void molecule::DetermineCenter(vector &Center)
 {
   atom *Walker = start;
@@ -788 +805 @@
   
   do {
-    CenterOfGravity.Zero();
+    Center.Zero();
     flag = true;
     while (Walker->next != end) {
@@ -815 +832 @@
         TestVector.AddVector(&TranslationVector);
         TestVector.MatrixMultiplication(matrix);
-        CenterOfGravity.AddVector(&TestVector);
+        Center.AddVector(&TestVector);
         cout << Verbose(1) << "Vector is: ";
         TestVector.Output((ofstream *)&cout);
@@ -828 +845 @@
             TestVector.AddVector(&TranslationVector);
             TestVector.MatrixMultiplication(matrix);
-            CenterOfGravity.AddVector(&TestVector);
+            Center.AddVector(&TestVector);
             cout << Verbose(1) << "Hydrogen Vector is: ";
             TestVector.Output((ofstream *)&cout);
@@ -838 +855 @@
     }
   } while (!flag);
-  Free((void **)&matrix, "molecule::DetermineCenterOfGravity: *matrix");
-  CenterOfGravity.Scale(1./(double)AtomCount);
+  Free((void **)&matrix, "molecule::DetermineCenter: *matrix");
+  Center.Scale(1./(double)AtomCount);
+};
+
+/** Transforms/Rotates the given molecule into its principal axis system.
+ * \param *out output stream for debugging
+ * \param DoRotate whether to rotate (true) or only to determine the PAS.
+ */
+void molecule::PrincipalAxisSystem(ofstream *out, bool DoRotate)
+{
+        atom *ptr = start;  // start at first in list
+        double InertiaTensor[NDIM*NDIM];
+        vector *CenterOfGravity = DetermineCenterOfGravity(out);
+
+        CenterGravity(out, CenterOfGravity);
+        
+        // reset inertia tensor 
+        for(int i=0;i<NDIM*NDIM;i++)
+                InertiaTensor[i] = 0.;
+        
+        // sum up inertia tensor
+        while (ptr->next != end) {
+                ptr = ptr->next;
+                vector x;
+                x.CopyVector(&ptr->x);
+                //x.SubtractVector(CenterOfGravity);
+                InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
+                InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
+                InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
+                InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
+                InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
+                InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
+                InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
+                InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
+                InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+        }
+        // print InertiaTensor for debugging
+        *out << "The inertia tensor is:" << endl;
+        for(int i=0;i<NDIM;i++) {
+          for(int j=0;j<NDIM;j++)
+            *out << InertiaTensor[i*NDIM+j] << " ";
+          *out << endl;
+        }
+        *out << endl;
+        
+        // diagonalize to determine principal axis system
+        gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
+        gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
+        gsl_vector *eval = gsl_vector_alloc(NDIM);
+        gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
+        gsl_eigen_symmv(&m.matrix, eval, evec, T);
+        gsl_eigen_symmv_free(T);
+        gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
+        
+        for(int i=0;i<NDIM;i++) {
+                *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
+                *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
+        }
+        
+        // check whether we rotate or not
+        if (DoRotate) {
+          *out << Verbose(1) << "Transforming molecule into PAS ... "; 
+          // the eigenvectors specify the transformation matrix
+          ptr = start;
+          while (ptr->next != end) {
+            ptr = ptr->next;
+            ptr->x.MatrixMultiplication(evec->data);
+          }
+          *out << "done." << endl;
+
+          // summing anew for debugging (resulting matrix has to be diagonal!)
+          // reset inertia tensor 
+    for(int i=0;i<NDIM*NDIM;i++)
+      InertiaTensor[i] = 0.;
+    
+    // sum up inertia tensor
+    ptr = start;
+    while (ptr->next != end) {
+      ptr = ptr->next;
+      vector x;
+      x.CopyVector(&ptr->x);
+      //x.SubtractVector(CenterOfGravity);
+      InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
+      InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
+      InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
+      InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
+      InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
+      InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
+      InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
+      InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
+      InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+    }
+    // print InertiaTensor for debugging
+    *out << "The inertia tensor is:" << endl;
+    for(int i=0;i<NDIM;i++) {
+      for(int j=0;j<NDIM;j++)
+        *out << InertiaTensor[i*NDIM+j] << " ";
+      *out << endl;
+    }
+    *out << endl;
+        }
+        
+        // free everything
+        delete(CenterOfGravity);
+        gsl_vector_free(eval);
+        gsl_matrix_free(evec);
 };
 
@@ -1234 +1355 @@
   if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
     *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
-    Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
+    Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
     while (Subgraphs->next != NULL) {
       Subgraphs = Subgraphs->next;
@@ -1248 +1369 @@
   }
   return No;
+};
+
+/** Determines the volume of a cluster.
+ * Determines first the convex envelope, then tesselates it and calculates its volume.
+ * \param *out output stream for debugging
+ * \param IsAngstroem for the units on output
+ */
+double molecule::VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem) 
+{
+        atom *Walker = NULL;
+        
+        // 1. calculate center of gravity
+        vector *CenterOfGravity = DetermineCenterOfGravity(out);
+        
+        // 2. calculate distance of each atom and sort into hash table
+//      map<double, int> DistanceSet;
+//      
+//      Walker = start;
+//      while (Walker->next != end) {
+//              Walker = Walker->next;
+//              DistanceSet.insert( pair<double, int>(ptr->x.distance(CenterOfGravity), ptr->nr) ); 
+//      }
+        
+        // 3. Find all points on the boundary
+        Boundaries BoundaryPoints[NDIM];        // first is alpha, second is (r, nr)
+        BoundariesTestPair BoundaryTestPair;
+        vector AxisVector;
+        // 3a. Go through every axis
+        for (int axis=0; axis<NDIM; axis++)  {
+                AxisVector.Zero();
+                AxisVector.x[axis] = 1.;
+                // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+                Walker = start;
+                while (Walker->next != end) {
+                        Walker = Walker->next;
+                        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) ) );
+                        if (BoundaryTestPair.second) { // successfully inserted
+                        } else { // same point exists, check first r, then distance of original vectors to centers of gravity
+                                *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
+                                *out << Verbose(2) << "Present vector: ";
+                                BoundaryTestPair.first->second.second->x.Output(out);
+                                *out << endl;
+                                *out << Verbose(2) << "New vector: ";
+                                Walker->x.Output(out);
+                                *out << endl;
+                                double tmp = ProjectedVector.Norm(); 
+                                if (tmp > BoundaryTestPair.first->second.first) {
+                                        BoundaryTestPair.first->second.first = tmp;
+                                        BoundaryTestPair.first->second.second = Walker; 
+                                        *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)) {
+                                                BoundaryTestPair.first->second.second = Walker;
+                                                *out << Verbose(2) << "Keeping new vector." << endl;
+                                        } else {
+                                                *out << Verbose(2) << "Keeping present vector." << endl;
+                                        }
+                                } else {
+                                                *out << Verbose(2) << "Keeping present vector." << 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
+                        flag = false;
+                        Boundaries::iterator left = BoundaryPoints[axis].end();
+                        Boundaries::iterator right = BoundaryPoints[axis].end();
+                        for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+                                // set neighbours correctly
+                                if (runner == BoundaryPoints[axis].begin()) {
+                                        left = BoundaryPoints[axis].end();
+                                }       else {
+                                        left = runner;
+                                }
+        left--;
+        right = runner;
+        right++;
+                                if (right == BoundaryPoints[axis].end()) {
+                                        right = BoundaryPoints[axis].begin();
+                                }
+                                // check distance
+                                if (runner->second.first < (left->second.first + right->second.first)/2. ) {
+                                        // throw out point
+                                        BoundaryPoints[axis].erase(runner);
+                                        flag = true;
+                                }
+                        }
+                } while (flag);
+        }       
+        // 3d. put into boundary set only those points appearing in each of the NDIM sets
+        int *AtomList = new int[AtomCount];
+        for (int j=0; j<AtomCount; j++) // reset list
+                AtomList[j] = 0;
+        for (int axis=0; axis<NDIM; axis++)  {  // increase whether it's present
+                for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+                        AtomList[runner->second.second->nr]++;
+                }
+        }
+        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)
+                                BoundaryPoints[axis].erase(runner);
+                }
+        }
+
+        *out << Verbose(2) << "I found " << BoundaryPointCount << " points on the convex boundary." << endl;
+        // now we have the whole set of edge points in the BoundaryList
+        
+        // 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;
+                        }
+                }
+        }
+        
+        // then create each line
+  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++;
+                        }
+      }
+                }
+        }
+        *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;
+                int nr;
+        } TriangleList[2 * 3 * AtomCount];      // each line can be part in at most 2 triangles
+        int TriangleCount = 0;
+        for (int i=0; i<LineCount; i++) {       // 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(2) << "I created " << TriangleCount << " triangles." << endl;
+
+        // 6. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+        double volume = 0.;
+        double PyramidVolume = 0.;
+        double G,h;
+        vector x,y;
+        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));
+                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)
+                volume += PyramidVolume; 
+        }
+        *out << Verbose(1) << "The summed volume is " << volume << " " << (IsAngstroem ? "A^3" : "a_0^3") << "." << endl;
+
+        // free reference lists
+        delete[](AtomList);
+  for(int i=0;i<AtomCount;i++)
+    delete[](AtomLines[i]);
+  delete[](AtomLines);
+
+        return volume;
 };
 
@@ -1490 +1861 @@
  * We use the algorithm from [Even, Graph Algorithms, p.62].
  * \param *out output stream for debugging
- * \param ReturnStack true - return pointer to atom stack of separable components, false - return NULL
  * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found
  * \return list of each disconnected subgraph as an individual molecule class structure
  */
-MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize)
+MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, int *&MinimumRingSize)
 {
   class StackClass<atom *> *AtomStack;
@@ -2006 +2376 @@
  * \param *path path to file
  * \param *FragmentList empty, filled on return
- * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
  */
-bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList, bool IsAngstroem)
+bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
 {
   bool status = true;
@@ -2186 +2555 @@
 {
   ifstream File;
-  stringstream line;
+  stringstream filename;
   bool status = true;
   char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
   
-  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
-  File.open(line.str().c_str(), ios::out);
+  filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  File.open(filename.str().c_str(), ios::out);
   *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
   if (File != NULL) {
@@ -2493 +2862 @@
 
   // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
-  Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
+  Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
   // fill the bond structure of the individually stored subgraphs
   Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
 
   // ===== 3. if structure still valid, parse key set file and others =====
-  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList, configuration->GetIsAngstroem());
+  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
 
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
@@ -3421 +3790 @@
   KeyStack AtomStack;
   atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
-  KeySet::iterator runner;
   int RootKeyNr = FragmentSearch.Root->nr;
   int Counter = FragmentSearch.FragmentCounter;
@@ -3816 +4184 @@
 {
   Graph ***FragmentLowerOrdersList = NULL;
-  int Order, NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
-  int counter = 0;
+  int NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
+  int counter = 0, Order;
   int UpgradeCount = RootStack.size();
   KeyStack FragmentRootStack;
@@ -3831 +4199 @@
 
   // initialise the fragments structure
-  FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
-  FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
   FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
   FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
@@ -4059 +4425 @@
   if (result) {
     *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
-    DetermineCenterOfGravity(CenterOfGravity);
-    OtherMolecule->DetermineCenterOfGravity(OtherCenterOfGravity);
+    DetermineCenter(CenterOfGravity);
+    OtherMolecule->DetermineCenter(OtherCenterOfGravity);
     *out << Verbose(5) << "Center of Gravity: ";
     CenterOfGravity.Output(out);

src/molecules.hpp

@@ -13 +13 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_matrix.h>
+#include <gsl/gsl_eigen.h>
 #include <gsl/gsl_heapsort.h>
 
@@ -43 +44 @@
 #define KeySetTestPair pair<KeySet::iterator, bool>
 #define GraphTestPair pair<Graph::iterator, bool>
+
+#define DistanceNrPair pair< double, atom* >
+#define Boundaries map <double, DistanceNrPair >
+#define BoundariesPair pair<double, DistanceNrPair >
+#define BoundariesTestPair pair<Boundaries::iterator, bool>
+#define LinetoAtomMap map < int, atom * > 
+#define LinetoAtomPair pair < int, atom * > 
+#define LinetoAtomTestPair pair < LinetoAtomMap::iterator, bool> 
 
 struct KeyCompare
@@ -265 +274 @@
   void CenterEdge(ofstream *out, vector *max); 
   void CenterOrigin(ofstream *out, vector *max); 
-  void CenterGravity(ofstream *out, vector *max); 
+  void CenterGravity(ofstream *out, vector *max);
   void Translate(const vector *x);
   void Mirror(const vector *x);
   void Align(vector *n);
   void Scale(double **factor);
-  void DetermineCenterOfGravity(vector &CenterOfGravity);
+  void DetermineCenter(vector &center);
+  vector * DetermineCenterOfGravity(ofstream *out);
   void SetBoxDimension(vector *dim);
   double * ReturnFullMatrixforSymmetric(double *cell_size);
   void ScanForPeriodicCorrection(ofstream *out);
-
+        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
+        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
+        
   bool CheckBounds(const vector *x) const;
   void GetAlignVector(struct lsq_params * par) const;
@@ -283 +295 @@
   
   // Graph analysis
-  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize);
+  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, int *&MinimumRingSize);
   void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
   bond * FindNextUnused(atom *vertex);
@@ -303 +315 @@
   bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
   bool StoreOrderAtSiteFile(ofstream *out, char *path);
-  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList, bool IsAngstroem);
+  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
   bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
   bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);