Changes in / [8a3e45:d250b2]

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• boundary.cpp (deleted)
• boundary.hpp (deleted)
• 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 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
+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
 
 bin_PROGRAMS = molecuilder joiner analyzer

src/builder.cpp

@@ -52 +52 @@
 #include "helpers.hpp"
 #include "molecules.hpp"
-#include "boundary.hpp"
 
 /********************************************** Submenu routine **************************************/
@@ -485 +484 @@
  * \param *mol the molecule with all the atoms
  */
-static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
+static void MeasureAtoms(periodentafel *periode, molecule *mol)
 {
   atom *first, *second, *third;
@@ -497 +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;
@@ -556 +553 @@
       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;
   }
 };
@@ -736 +720 @@
   int ExitFlag = 0;
   int j;
-  double volume = 0.;
   enum ConfigStatus config_present = absent;
   clock_t start,end;
@@ -759 +742 @@
             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;
@@ -790 +768 @@
             argptr+=1;
             break;
-          case 'n':
-            cout << "I won't parse trajectories." << endl;
-            configuration.FastParsing = true;
           default:   // no match? Step on
             argptr++;
@@ -805 +780 @@
       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
@@ -972 +945 @@
               }
               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
-              if (argv[argptr][0] != '-') // if it started with a '-' we've already made a step!
-                argptr++;
+              argptr++;
               break;
           }
@@ -1271 +1151 @@
 
       case 'l': // measure distances or angles
-        MeasureAtoms(periode, mol, &configuration);
+        MeasureAtoms(periode, mol);
         break;
 
@@ -1284 +1164 @@
         mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
         //mol->CreateListOfBondsPerAtom((ofstream *)&cout);
-        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, MinimumRingSize);
+        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
         while (Subgraphs->next != NULL) {
           Subgraphs = Subgraphs->next;
@@ -1340 +1220 @@
   
   // save element data base
-  if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
+  if (periode->StorePeriodentafel()) //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;
@@ -456 +455 @@
   string zeile;
   string dummy;
-  element *elementhash[MAX_ELEMENTS];
-  char name[MAX_ELEMENTS];
-  char keyword[MAX_ELEMENTS];
-  int Z, No[MAX_ELEMENTS];
+  element *elementhash[128];
+  char name[128];
+  char keyword[128];
+  int Z, No;
   int verbose = 0;
   
@@ -633 +632 @@
   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[i], 1, critical);
+    ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No, 1, critical);
     ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
     elementhash[i] = periode->FindElement(Z);
-    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++) {
+    for(int j=0;j<No;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
-      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;
-      }
+      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);
@@ -665 +656 @@
       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, (int)FastParsing, 7, 1, double_type, &neues->v.x[2], repetition,optional))
+      if(!ParseForParameter(verbose,file, keyword, 0, 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"    //!< 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 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 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", PathBackup, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    //strcpy(PathBackup, configuration->configpath);
-    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
+    strcpy(PathBackup, configuration->configpath);
+    sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     
     // center on edge
@@ -251 +251 @@
     for (int k=0;k<NDIM;k++) {
       j += k+1;
-      BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
+      BoxDimension.x[k] = 5.;
       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
@@ -252 +253 @@
     case 2:
       // determine two other bonds (warning if there are more than two other) plus valence sanity check
-      for (int i=0;i<NumBond;i++) {
+      for (i=0;i<NumBond;i++) {
         if (BondList[i] != TopBond) {
           if (FirstBond == NULL) {
@@ -311 +312 @@
       FirstOtherAtom->x.Zero();
       SecondOtherAtom->x.Zero();
-      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
+      for(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));
@@ -318 +319 @@
       SecondOtherAtom->x.Scale(&BondRescale);
       //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
-      for(int i=NDIM;i--;) { // and make relative to origin atom
+      for(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];
@@ -438 +439 @@
   int NumberOfAtoms = 0; // atom number in xyz read
   int i, j; // loop variables
-  atom *Walker = NULL;  // pointer to added atom
+  atom *first = NULL;  // pointer to added atom
   char shorthand[3];  // shorthand for atom name
   ifstream xyzfile;   // xyz file
@@ -456 +457 @@
 
   for(i=0;i<NumberOfAtoms;i++){
-    Walker = new atom;
+    first = new atom;
     getline(xyzfile,line,'\n');
     istringstream *item = new istringstream(line);
@@ -465 +466 @@
     *item >> x[1];
     *item >> x[2];
-    Walker->type = elemente->FindElement(shorthand);
-    if (Walker->type == NULL) {
+    first->type = elemente->FindElement(shorthand);
+    if (first->type == NULL) {
       cerr << "Could not parse the element at line: '" << line << "', setting to H.";
-      Walker->type = elemente->FindElement(1);
+      first->type = elemente->FindElement(1);
     }
     for(j=NDIM;j--;)
-      Walker->x.x[j] = x[j];
-     AddAtom(Walker);  // add to molecule
+      first->x.x[j] = x[j];
+     AddAtom(first);  // add to molecule
     delete(item);
   }
@@ -709 +710 @@
 }; 
 
-/** Returns vector pointing to center of gravity.
+/** Centers the center of gravity of the atoms at (0,0,0).
  * \param *out output stream for debugging
- * \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;
+ * \param *center return vector for translation vector
+ */
+void molecule::CenterGravity(ofstream *out, vector *center)
+{
   double Num = 0;
-        
-        a->Zero();
-
+  atom *ptr = start->next;  // start at first in list
+  vector tmp;
+  
+  for(int i=NDIM;i--;) // zero center vector
+    center->x[i] = 0.;
+    
   if (ptr != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
@@ -728 +729 @@
       tmp.CopyVector(&ptr->x);
       tmp.Scale(ptr->type->mass);  // scale by mass
-      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 {
+      center->AddVector(&tmp);       
+    }
+    center->Scale(-1./Num); // divide through total mass (and sign for direction)
     Translate(center);
   }
@@ -792 +775 @@
 };
 
-/** Determines center of molecule (yet not considering atom masses).
- * \param Center reference to return vector
- */
-void molecule::DetermineCenter(vector &Center)
+/** Determines center of gravity (yet not considering atom masses).
+ * \param CenterOfGravity reference to return vector
+ */
+void molecule::DetermineCenterOfGravity(vector &CenterOfGravity)
 {
   atom *Walker = start;
@@ -805 +788 @@
   
   do {
-    Center.Zero();
+    CenterOfGravity.Zero();
     flag = true;
     while (Walker->next != end) {
@@ -832 +815 @@
         TestVector.AddVector(&TranslationVector);
         TestVector.MatrixMultiplication(matrix);
-        Center.AddVector(&TestVector);
+        CenterOfGravity.AddVector(&TestVector);
         cout << Verbose(1) << "Vector is: ";
         TestVector.Output((ofstream *)&cout);
@@ -845 +828 @@
             TestVector.AddVector(&TranslationVector);
             TestVector.MatrixMultiplication(matrix);
-            Center.AddVector(&TestVector);
+            CenterOfGravity.AddVector(&TestVector);
             cout << Verbose(1) << "Hydrogen Vector is: ";
             TestVector.Output((ofstream *)&cout);
@@ -855 +838 @@
     }
   } while (!flag);
-  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);
+  Free((void **)&matrix, "molecule::DetermineCenterOfGravity: *matrix");
+  CenterOfGravity.Scale(1./(double)AtomCount);
 };
 
@@ -1355 +1234 @@
   if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
     *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
-    Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
+    Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
     while (Subgraphs->next != NULL) {
       Subgraphs = Subgraphs->next;
@@ -1370 +1249 @@
   return No;
 };
+
 /** Returns Shading as a char string.
  * \param color the Shading
@@ -1573 +1453 @@
             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
@@ -1579 +1459 @@
               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++;
@@ -1591 +1471 @@
           *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; 
@@ -1610 +1490 @@
  * 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, int *&MinimumRingSize)
+MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize)
 {
   class StackClass<atom *> *AtomStack;
@@ -2125 +2006 @@
  * \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 molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList, bool IsAngstroem)
 {
   bool status = true;
@@ -2304 +2186 @@
 {
   ifstream File;
-  stringstream filename;
+  stringstream line;
   bool status = true;
   char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
   
-  filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
-  File.open(filename.str().c_str(), ios::out);
+  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
+  File.open(line.str().c_str(), ios::out);
   *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
   if (File != NULL) {
@@ -2611 +2493 @@
 
   // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
-  Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
+  Subgraphs = DepthFirstSearchAnalysis(out, false, 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);
+  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList, configuration->GetIsAngstroem());
 
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
@@ -3131 +3013 @@
 //          *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
@@ -3539 +3421 @@
   KeyStack AtomStack;
   atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
+  KeySet::iterator runner;
   int RootKeyNr = FragmentSearch.Root->nr;
   int Counter = FragmentSearch.FragmentCounter;
@@ -3933 +3816 @@
 {
   Graph ***FragmentLowerOrdersList = NULL;
-  int NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
-  int counter = 0, Order;
+  int Order, NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
+  int counter = 0;
   int UpgradeCount = RootStack.size();
   KeyStack FragmentRootStack;
@@ -3948 +3831 @@
 
   // 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");
@@ -4174 +4059 @@
   if (result) {
     *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
-    DetermineCenter(CenterOfGravity);
-    OtherMolecule->DetermineCenter(OtherCenterOfGravity);
+    DetermineCenterOfGravity(CenterOfGravity);
+    OtherMolecule->DetermineCenterOfGravity(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>
 
@@ -40 +39 @@
 #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>
@@ -266 +265 @@
   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 DetermineCenter(vector &center);
-  vector * DetermineCenterOfGravity(ofstream *out);
+  void DetermineCenterOfGravity(vector &CenterOfGravity);
   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;
@@ -287 +283 @@
   
   // Graph analysis
-  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, int *&MinimumRingSize);
+  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize);
   void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
   bond * FindNextUnused(atom *vertex);
@@ -307 +303 @@
   bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
   bool StoreOrderAtSiteFile(ofstream *out, char *path);
-  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
+  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList, bool IsAngstroem);
   bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
   bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
@@ -396 +392 @@
     char *configpath;
     char *configname;
-    bool FastParsing;
     
     private:

src/periodentafel.cpp

@@ -217 +217 @@
 
   // fill valence DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, path, MAXSTRINGSIZE);
+  strncpy(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -227 +226 @@
         infile >> FindElement((int)tmp)->Valence;
         infile >> ws;
-        //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->Valence << " valence electrons." << endl;
+        //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->Valence << " valence electrons." << endl;
     }
     infile.close();
@@ -235 +234 @@
 
   // fill valence DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, path, MAXSTRINGSIZE);
+  strncpy(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -245 +243 @@
       infile >> FindElement((int)tmp)->NoValenceOrbitals;
       infile >> ws;
-      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
     }
     infile.close();
@@ -253 +251 @@
   
   // fill H-BondDistance DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, path, MAXSTRINGSIZE);
+  strncpy(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -266 +263 @@
       infile >> ptr->HBondDistance[2];
         infile >> ws;
-      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << find_element((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
     }
     infile.close();
@@ -274 +271 @@
   
   // fill H-BondAngle DB per element
-  strncpy(filename, path, MAXSTRINGSIZE);
-  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
-  strncat(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, path, MAXSTRINGSIZE);
+  strncpy(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -293 +289 @@
     otherstatus = false;
   
-  if (!otherstatus)
-    cerr << "ERROR: Something went wrong while parsing the databases!" << endl;
-  
   return status;
 };
@@ -305 +298 @@
   bool result = true;
   ofstream f;
-  char file[MAXSTRINGSIZE];
-  
-  if (filename == STANDARDELEMENTSDB)
-    f.open(file);
+  
+  if (filename == NULL)
+    f.open("elements.db");
   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.
@@ -115 +111 @@
 };
 
-/** 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
@@ -133 +117 @@
 double vector::Projection(const vector *y) const
 {
-  return (ScalarProduct(y));
+  return (ScalarProduct(y)/Norm()/y->Norm());
 };
 
@@ -166 +150 @@
 };
 
-/** 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$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
+ * \return \f$\frac{\langle x, y \rangle}{|x||y|}\f$
  */
 double vector::Angle(vector *y) const
 {
-  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
+  return (this->ScalarProduct(y)/(this->Norm()*y->Norm()));
 };
 
@@ -387 +354 @@
 {
   double projection;
-  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
+  projection = ScalarProduct(n)/((vector *)n)->ScalarProduct(n);    // remove constancy from n (keep as logical one)
   // withdraw projected vector twice from original one
   cout << Verbose(1) << "Vector: ";
@@ -418 +385 @@
     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]);
@@ -452 +419 @@
     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]);
@@ -491 +458 @@
  * \return true - success, false - failure (null vector given)
  */
-bool vector::GetOneNormalVector(const vector *GivenVector)
+bool vector::GetOneNormalVector(const vector *vector)
 {
   int Components[NDIM]; // contains indices of non-zero components
@@ -499 +466 @@
 
   cout << Verbose(4);
-  GivenVector->Output((ofstream *)&cout);
+  vector->Output((ofstream *)&cout);
   cout << endl;
   for (j=NDIM;j--;)
@@ -505 +472 @@
   // find two components != 0
   for (j=0;j<NDIM;j++)
-    if (fabs(GivenVector->x[j]) > MYEPSILON)
+    if (fabs(vector->x[j]) > MYEPSILON)
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
@@ -512 +479 @@
     case 3:  // threecomponent system
     case 2:  // two component system
-      norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
+      norm = sqrt(1./(vector->x[Components[1]]*vector->x[Components[1]]) + 1./(vector->x[Components[0]]*vector->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./GivenVector->x[Components[1]] / norm;
-      x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
+      x[Components[1]] = -1./vector->x[Components[1]] / norm;
+      x[Components[0]] = 1./vector->x[Components[0]] / norm;
       return true;
       break;
@@ -531 +498 @@
 };
 
-/** 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
@@ -568 +519 @@
      gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *s = NULL;
-   gsl_vector *ss, *y;
+   gsl_vector *ss, *x;
    gsl_multimin_function minex_func;
  
@@ -577 +528 @@
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
-   y = gsl_vector_alloc (np);
+   x = gsl_vector_alloc (np);
  
    /* Set all step sizes to 1 */
@@ -587 +538 @@
         
          for (i=NDIM;i--;)
-                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
+                gsl_vector_set(x, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 
          
    /* Initialize method and iterate */
@@ -595 +546 @@
  
    s = gsl_multimin_fminimizer_alloc (T, np);
-   gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
+   gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
  
    do
@@ -624 +575 @@
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
-   gsl_vector_free(y);
+   gsl_vector_free(x);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free (s);

src/vector.hpp

@@ -10 +10 @@
 
   vector();
-  vector(double x1, double x2, double x3);
   ~vector();
 
@@ -24 +23 @@
   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);
@@ -39 +35 @@
   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);