Changes in / [d250b2:8a3e45]

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• boundary.cpp (added)
• boundary.hpp (added)
• builder.cpp (modified) (12 diffs)
• config.cpp (modified) (4 diffs)
• defs.hpp (modified) (1 diff)
• helpers.hpp (modified) (1 diff)
• moleculelist.cpp (modified) (2 diffs)
• molecules.cpp (modified) (28 diffs)
• molecules.hpp (modified) (6 diffs)
• periodentafel.cpp (modified) (9 diffs)
• vector.cpp (modified) (17 diffs)
• vector.hpp (modified) (3 diffs)

src/Makefile.am ¶

@@ -1 @@
-SOURCE = atom.cpp bond.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
-HEADER = defs.hpp helpers.hpp molecules.hpp stackclass.hpp vector.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
+HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp stackclass.hpp vector.hpp
 
 bin_PROGRAMS = molecuilder joiner analyzer

src/builder.cpp ¶

@@ -52 +52 @@
 #include "helpers.hpp"
 #include "molecules.hpp"
+#include "boundary.hpp"
 
 /********************************************** Submenu routine **************************************/
@@ -484 +485 @@
  * \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 +497 @@
   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 +556 @@
       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.";
+        VolumeOfConvexEnvelope((ofstream *)&cout, configuration, NULL, mol);
+        break;
   }
 };
@@ -720 +736 @@
   int ExitFlag = 0;
   int j;
+  double volume = 0.;
   enum ConfigStatus config_present = absent;
   clock_t start,end;
@@ -742 +759 @@
             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;
+            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
             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-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-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;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << 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-v/-V\t\tGives version information." << endl;
             cout << "Note: config files must not begin with '-' !" << endl;
@@ -768 +790 @@
             argptr+=1;
             break;
+          case 'n':
+            cout << "I won't parse trajectories." << endl;
+            configuration.FastParsing = true;
           default:   // no match? Step on
             argptr++;
@@ -780 +805 @@
       cout << Verbose(0) << "Element list loaded successfully." << endl;
       periode->Output((ofstream *)&cout);
-    } else
+    } else {
       cout << Verbose(0) << "Element list loading failed." << endl;
+      return 1;
+    }
     
     // 3. Find config file name and parse if possible
@@ -945 +972 @@
               }
               break;
+            case 'm':
+              ExitFlag = 1;
+              j = atoi(argv[argptr++]);
+              if ((j<0) || (j>1)) {
+                cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
+                j = 0;
+              }
+              if (j)
+                cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
+              else
+                cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+              mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
+              break;
+            case 'o':
+              ExitFlag = 1;
+              cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+              VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
+              break;
+            case 'U':
+              volume = atof(argv[argptr++]);
+              cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;  
+            case 'u':
+              {
+                double density;
+                ExitFlag = 1;
+                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
+                density = atof(argv[argptr++]);
+                if (density < 1.0) {
+                  cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
+                  density = 1.3;
+                }
+//                for(int i=0;i<NDIM;i++) {
+//                  repetition[i] = atoi(argv[argptr++]);
+//                  if (repetition[i] < 1)
+//                    cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
+//                  repetition[i] = 1;
+//                }
+                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);
+              }
+              break;
+            case 'd':
+              for (int axis = 1; axis <= NDIM; axis++) {
+                int faktor = atoi(argv[argptr++]);
+                int count;
+                element ** Elements;
+                vector ** Vectors;
+                if (faktor < 1) {
+                  cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
+                  faktor = 1;
+                }
+                mol->CountAtoms((ofstream *)&cout);  // recount atoms
+                if (mol->AtomCount != 0) {  // if there is more than none
+                  count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
+                  Elements = new element *[count];
+                  Vectors = new vector *[count];
+                  j = 0;
+                  first = mol->start;
+                  while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
+                    first = first->next;
+                    Elements[j] = first->type;
+                    Vectors[j] = &first->x;
+                    j++;
+                  }
+                  if (count != j)
+                    cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+                  x.Zero();
+                  y.Zero();
+                  y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+                  for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
+                    x.AddVector(&y); // per factor one cell width further
+                    for (int k=count;k--;) { // go through every atom of the original cell
+                      first = new atom(); // create a new body
+                      first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
+                      first->x.AddVector(&x);      // translate the coordinates
+                      first->type = Elements[k];  // insert original element
+                      mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
+                    }
+                  }
+                  // free memory
+                  delete[](Elements);
+                  delete[](Vectors);
+                  // correct cell size
+                  if (axis < 0) { // if sign was negative, we have to translate everything
+                    x.Zero();
+                    x.AddVector(&y);
+                    x.Scale(-(faktor-1));
+                    mol->Translate(&x);
+                  }
+                  mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+                }
+              }
+              break;
             default:   // no match? Step on
-              argptr++;
+              if (argv[argptr][0] != '-') // if it started with a '-' we've already made a step!
+                argptr++;
               break;
           }
@@ -1151 +1271 @@
 
       case 'l': // measure distances or angles
-        MeasureAtoms(periode, mol);
+        MeasureAtoms(periode, mol, &configuration);
         break;
 
@@ -1164 +1284 @@
         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;
@@ -1220 +1340 @@
   
   // save element data base
-  if (periode->StorePeriodentafel()) //ElementsFileName
+  if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
     cout << Verbose(0) << "Saving of elements.db successful." << endl;
   else

src/config.cpp ¶

@@ -24 +24 @@
   configname[0]='\0';
   
+  FastParsing = false;
   ProcPEGamma=8;
   ProcPEPsi=1;
@@ -455 +456 @@
   string zeile;
   string dummy;
-  element *elementhash[128];
-  char name[128];
-  char keyword[128];
-  int Z, No;
+  element *elementhash[MAX_ELEMENTS];
+  char name[MAX_ELEMENTS];
+  char keyword[MAX_ELEMENTS];
+  int Z, No[MAX_ELEMENTS];
   int verbose = 0;
   
@@ -632 +633 @@
   if (!ParseForParameter(verbose,file,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
     config::StructOpt = 0;
+  // prescan number of ions per type
+  cout << Verbose(0) << "Prescanning ions per type: " << endl;
   for (int i=0; i < config::MaxTypes; i++) {
     sprintf(name,"Ion_Type%i",i+1);
-    ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No, 1, critical);
+    ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
     ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
     elementhash[i] = periode->FindElement(Z);
-    for(int j=0;j<No;j++) {
+    cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl; 
+  }
+  for (int i=0; i < config::MaxTypes; i++) {
+    sprintf(name,"Ion_Type%i",i+1);
+    for(int j=0;j<No[i];j++) {
       int repetition = 1; // which repeated keyword shall be read
       sprintf(keyword,"%s_%i",name, j+1);
       atom *neues = new atom();
       // then parse for each atom the coordinates as often as present
-      while (ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition, (repetition == 1) ? critical : optional) &&
-        ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition, (repetition == 1) ? critical : optional) &&
-        ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition, (repetition == 1) ? critical : optional))
-        repetition++;
-      repetition--;
-      //cout << "Found " << repetition << " times of the keyword " << keyword << "." << endl;
+      if (!FastParsing) {
+        while (ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition, (repetition == 1) ? critical : optional) &&
+          ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition, (repetition == 1) ? critical : optional) &&
+          ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition, (repetition == 1) ? critical : optional))
+          repetition++;
+        repetition--;
+        //cout << "Found " << repetition << " times of the keyword " << keyword << "." << endl;
+      }
       ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
       ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
@@ -656 +665 @@
       if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
         neues->v.x[1] = 0.;
-      if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
+      if(!ParseForParameter(verbose,file, keyword, (int)FastParsing, 7, 1, double_type, &neues->v.x[2], repetition,optional))
         neues->v.x[2] = 0.;
       neues->type = elementhash[i]; // find element type

src/defs.hpp ¶

@@ -43 +43 @@
 
 // various standard filenames
-#define DEFAULTCONFIG "main_pcp_linux"
-#define KEYSETFILE "KeySets.dat"
-#define ADJACENCYFILE "Adjacency.dat"
-#define TEFACTORSFILE "TE-Factors.dat"
-#define FORCESFILE "Forces-Factors.dat"
-#define ORDERATSITEFILE "OrderAtSite.dat"
-#define ENERGYPERFRAGMENT "EnergyPerFragment"
-#define FRAGMENTPREFIX "BondFragment"
-#define STANDARDCONFIG "unknown.conf"
-#define STANDARDELEMENTSDB "elements.db"
-#define STANDARDVALENCEDB "valence.db"
-#define STANDARDORBITALDB "orbitals.db"
-#define STANDARDHBONDDISTANCEDB "Hbonddistance.db"
-#define STANDARDHBONDANGLEDB "Hbondangle.db"
+#define DEFAULTCONFIG "main_pcp_linux"    //!< default filename of config file
+#define CONVEXENVELOPE "ConvexEnvelope.dat"    //!< default filename of convex envelope tecplot data file
+#define KEYSETFILE "KeySets.dat"    //!< default filename of BOSSANOVA key sets file
+#define ADJACENCYFILE "Adjacency.dat"    //!< default filename of BOSSANOVA adjacancy file
+#define TEFACTORSFILE "TE-Factors.dat"    //!< default filename of BOSSANOVA total energy factors file
+#define FORCESFILE "Forces-Factors.dat"    //!< default filename of BOSSANOVA force factors file
+#define ORDERATSITEFILE "OrderAtSite.dat"    //!< default filename of BOSSANOVA Bond Order at each atom file
+#define ENERGYPERFRAGMENT "EnergyPerFragment"    //!< default filename of BOSSANOVA Energy contribution Per Fragment file
+#define FRAGMENTPREFIX "BondFragment"    //!< default filename prefix of BOSSANOVA fragment config and directories 
+#define STANDARDCONFIG "unknown.conf"    //!< default filename of standard config file
+#define STANDARDELEMENTSDB "elements.db"    //!< default filename of elements data base with masses, Z, VanDerWaals radii, ...
+#define STANDARDVALENCEDB "valence.db"    //!< default filename of valence number per element database
+#define STANDARDORBITALDB "orbitals.db"    //!< default filename of orbitals per element database
+#define STANDARDHBONDDISTANCEDB "Hbonddistance.db"    //!< default filename of typial bond distance to hydrogen database
+#define STANDARDHBONDANGLEDB "Hbondangle.db"    //!< default filename of typial bond angle to hydrogen database
+
+// some values
+#define SOLVENTDENSITY_A 0.6022142
+#define SOLVENTDENSITY_a0 0.089238936
+
 
 #define UPDATECOUNT 10  //!< update ten sites per BOSSANOVA interval

src/helpers.hpp ¶

@@ -118 +118 @@
 };
 
-
-
 /******************************** Some templates for list management ***********************************/
 

src/moleculelist.cpp ¶

@@ -240 +240 @@
     *out << endl;
     // prepare output of config file
-    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.conf", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    strcpy(PathBackup, configuration->configpath);
-    sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    //strcpy(PathBackup, configuration->configpath);
+    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     
     // center on edge
@@ -251 +251 @@
     for (int k=0;k<NDIM;k++) {
       j += k+1;
-      BoxDimension.x[k] = 5.;
+      BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
       ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
     }

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;
@@ -1249 +1370 @@
   return No;
 };
-
 /** Returns Shading as a char string.
  * \param color the Shading
@@ -1453 +1573 @@
             for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
               NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-                  //*out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+                  *out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
                   if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check NoBonds of other atom
                     // count valence of second partner
@@ -1459 +1579 @@
               for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
                 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
-              //*out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+              *out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
                     if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) // increase bond degree by one
                 ListOfBondsPerAtom[Walker->nr][i]->BondDegree++;
@@ -1471 +1591 @@
           *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << "." << endl;
                 
-//        // output bonds for debugging (if bond chain list was correctly installed)
-//        *out << Verbose(1) << endl << "From contents of bond chain list:";
-//        bond *Binder = first;
-//    while(Binder->next != last) {
-//      Binder = Binder->next;
-//                      *out << *Binder << "\t" << endl;
-//    }
-//    *out << endl;
+          // output bonds for debugging (if bond chain list was correctly installed)
+          *out << Verbose(1) << endl << "From contents of bond chain list:";
+          bond *Binder = first;
+    while(Binder->next != last) {
+      Binder = Binder->next;
+                        *out << *Binder << "\t" << endl;
+    }
+    *out << endl;
   } else
         *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl; 
@@ -1490 +1610 @@
  * 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 +2125 @@
  * \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 +2304 @@
 {
   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 +2611 @@
 
   // ===== 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)
@@ -3013 +3131 @@
 //          *out << ", who has no son in this fragment molecule." << endl; 
 #ifdef ADDHYDROGEN
-//          *out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
+          //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
           Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem);
 #endif
@@ -3421 +3539 @@
   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 +3933 @@
 {
   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 +3948 @@
 
   // 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 +4174 @@
   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>
 
@@ -39 +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>
@@ -265 +266 @@
   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 +287 @@
   
   // 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 +307 @@
   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);
@@ -392 +396 @@
     char *configpath;
     char *configname;
+    bool FastParsing;
     
     private:

src/periodentafel.cpp ¶

@@ -217 +217 @@
 
   // fill valence DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -226 +227 @@
         infile >> FindElement((int)tmp)->Valence;
         infile >> ws;
-        //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->Valence << " valence electrons." << endl;
+        //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->Valence << " valence electrons." << endl;
     }
     infile.close();
@@ -234 +235 @@
 
   // fill valence DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -243 +245 @@
       infile >> FindElement((int)tmp)->NoValenceOrbitals;
       infile >> ws;
-      //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
     }
     infile.close();
@@ -251 +253 @@
   
   // fill H-BondDistance DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -263 +266 @@
       infile >> ptr->HBondDistance[2];
         infile >> ws;
-      //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
     }
     infile.close();
@@ -271 +274 @@
   
   // fill H-BondAngle DB per element
-  strncat(filename, path, MAXSTRINGSIZE);
-  strncpy(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -289 +293 @@
     otherstatus = false;
   
+  if (!otherstatus)
+    cerr << "ERROR: Something went wrong while parsing the databases!" << endl;
+  
   return status;
 };
@@ -298 +305 @@
   bool result = true;
   ofstream f;
-  
-  if (filename == NULL)
-    f.open("elements.db");
+  char file[MAXSTRINGSIZE];
+  
+  if (filename == STANDARDELEMENTSDB)
+    f.open(file);
   else
     f.open(filename);

src/vector.cpp ¶

@@ -13 +13 @@
  */
 vector::vector() { x[0] = x[1] = x[2] = 0.; };
+
+/** Constructor of class vector.
+ */
+vector::vector(double x1, double x2, double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
 
 /** Desctructor of class vector.
@@ -111 +115 @@
 };
 
+/** projects this vector onto plane defined by \a *y.
+ * \param *y array to normal vector of plane
+ * \return \f$\langle x, y \rangle\f$
+ */
+void vector::ProjectOntoPlane(const vector *y)
+{
+  vector tmp;
+  tmp.CopyVector(y);
+  tmp.Scale(Projection(y));
+  this->SubtractVector(&tmp);
+};
+
 /** Calculates the projection of a vector onto another \a *y.
  * \param *y array to second vector
@@ -117 +133 @@
 double vector::Projection(const vector *y) const
 {
-  return (ScalarProduct(y)/Norm()/y->Norm());
+  return (ScalarProduct(y));
 };
 
@@ -150 +166 @@
 };
 
+/** Zeros all components of this vector.
+ */
+void vector::One(double one)
+{
+  for (int i=NDIM;i--;)
+    this->x[i] = one;
+};
+
+/** Initialises all components of this vector.
+ */
+void vector::Init(double x1, double x2, double x3)
+{ 
+  x[0] = x1;
+  x[1] = x2;
+  x[2] = x3;
+};
+
 /** Calculates the angle between this and another vector.
  * \param *y array to second vector
- * \return \f$\frac{\langle x, y \rangle}{|x||y|}\f$
+ * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
  */
 double vector::Angle(vector *y) const
 {
-  return (this->ScalarProduct(y)/(this->Norm()*y->Norm()));
+  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
 };
 
@@ -354 +387 @@
 {
   double projection;
-  projection = ScalarProduct(n)/((vector *)n)->ScalarProduct(n);    // remove constancy from n (keep as logical one)
+  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
   // withdraw projected vector twice from original one
   cout << Verbose(1) << "Vector: ";
@@ -385 +418 @@
     return false;
   }
-  cout << Verbose(4) << "relative, first plane coordinates:";
-  x1.Output((ofstream *)&cout);
-  cout << endl;
-  cout << Verbose(4) << "second plane coordinates:";
-  x2.Output((ofstream *)&cout);
-  cout << endl;
+//  cout << Verbose(4) << "relative, first plane coordinates:";
+//  x1.Output((ofstream *)&cout);
+//  cout << endl;
+//  cout << Verbose(4) << "second plane coordinates:";
+//  x2.Output((ofstream *)&cout);
+//  cout << endl;
 
   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
@@ -419 +452 @@
     return false;
   }
-  cout << Verbose(4) << "relative, first plane coordinates:";
-  x1.Output((ofstream *)&cout);
-  cout << endl;
-  cout << Verbose(4) << "second plane coordinates:";
-  x2.Output((ofstream *)&cout);
-  cout << endl;
+//  cout << Verbose(4) << "relative, first plane coordinates:";
+//  x1.Output((ofstream *)&cout);
+//  cout << endl;
+//  cout << Verbose(4) << "second plane coordinates:";
+//  x2.Output((ofstream *)&cout);
+//  cout << endl;
 
   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
@@ -458 +491 @@
  * \return true - success, false - failure (null vector given)
  */
-bool vector::GetOneNormalVector(const vector *vector)
+bool vector::GetOneNormalVector(const vector *GivenVector)
 {
   int Components[NDIM]; // contains indices of non-zero components
@@ -466 +499 @@
 
   cout << Verbose(4);
-  vector->Output((ofstream *)&cout);
+  GivenVector->Output((ofstream *)&cout);
   cout << endl;
   for (j=NDIM;j--;)
@@ -472 +505 @@
   // find two components != 0
   for (j=0;j<NDIM;j++)
-    if (fabs(vector->x[j]) > MYEPSILON)
+    if (fabs(GivenVector->x[j]) > MYEPSILON)
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
@@ -479 +512 @@
     case 3:  // threecomponent system
     case 2:  // two component system
-      norm = sqrt(1./(vector->x[Components[1]]*vector->x[Components[1]]) + 1./(vector->x[Components[0]]*vector->x[Components[0]]));
+      norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
       x[Components[2]] = 0.;
       // in skp both remaining parts shall become zero but with opposite sign and third is zero
-      x[Components[1]] = -1./vector->x[Components[1]] / norm;
-      x[Components[0]] = 1./vector->x[Components[0]] / norm;
+      x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
+      x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
       return true;
       break;
@@ -498 +531 @@
 };
 
+/** Determines paramter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
+ * \param *A first plane vector
+ * \param *B second plane vector
+ * \param *C third plane vector
+ * \return scaling parameter for this vector
+ */
+double vector::CutsPlaneAt(vector *A, vector *B, vector *C)
+{
+//  cout << Verbose(3) << "For comparison: ";
+//  cout << "A " << A->Projection(this) << "\t"; 
+//  cout << "B " << B->Projection(this) << "\t"; 
+//  cout << "C " << C->Projection(this) << "\t"; 
+//  cout << endl;
+  return A->Projection(this);
+};
+
 /** Creates a new vector as the one with least square distance to a given set of \a vectors.
  * \param *vectors set of vectors
@@ -519 +568 @@
      gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *s = NULL;
-   gsl_vector *ss, *x;
+   gsl_vector *ss, *y;
    gsl_multimin_function minex_func;
  
@@ -528 +577 @@
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
-   x = gsl_vector_alloc (np);
+   y = gsl_vector_alloc (np);
  
    /* Set all step sizes to 1 */
@@ -538 +587 @@
         
          for (i=NDIM;i--;)
-                gsl_vector_set(x, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
+                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
          
    /* Initialize method and iterate */
@@ -546 +595 @@
  
    s = gsl_multimin_fminimizer_alloc (T, np);
-   gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
+   gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
  
    do
@@ -575 +624 @@
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
-   gsl_vector_free(x);
+   gsl_vector_free(y);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free (s);

src/vector.hpp ¶

@@ -10 +10 @@
 
   vector();
+  vector(double x1, double x2, double x3);
   ~vector();
 
@@ -23 +24 @@
   void CopyVector(const vector *y);
   void RotateVector(const vector *y, const double alpha);
+  void ProjectOntoPlane(const vector *y);
   void Zero(); 
+  void One(double one);
+  void Init(double x1, double x2, double x3);
   void Normalize();
   void Translate(const vector *x);
@@ -35 +39 @@
   void LinearCombinationOfVectors(const vector *x1, const vector *x2, const vector *x3, double *factors);
   
+  double CutsPlaneAt(vector *A, vector *B, vector *C);
   bool GetOneNormalVector(const vector *x1);
   bool MakeNormalVector(const vector *y1);