Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 14de4697703b4c1283d68496e457f6c1c3c794aa)
+++ src/builder.cpp	(revision 14de4697703b4c1283d68496e457f6c1c3c794aa)
@@ -0,0 +1,1170 @@
+/** \file builder.cpp
+ * 
+ * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
+ * The output is the complete configuration file for PCP for direct use.
+ * Features:
+ * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
+ * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
+ *    
+ */
+
+/*! \mainpage Molecuilder - a molecular set builder
+ * 
+ * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
+ * 
+ * \section about About the Program
+ * 
+ *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
+ *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
+ *  already constructed atoms.
+ * 
+ *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
+ *  molecular dynamics implementation.
+ * 
+ * \section install Installation
+ * 
+ *  Installation should without problems succeed as follows:
+ *  -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
+ *  -# make
+ *  -# make install
+ * 
+ *  Further useful commands are
+ *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
+ *  -# make doxygen-doc: Creates these html pages out of the documented source 
+ * 
+ * \section run Running
+ * 
+ *  The program can be executed by running: ./molecuilder
+ * 
+ *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
+ *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
+ *  later re-execution. 
+ * 
+ * \section ref References
+ *  
+ *  For the special configuration file format, see the documentation of pcp.
+ * 
+ */
+
+
+using namespace std;
+
+#include "helpers.hpp"
+#include "molecules.hpp"
+
+/********************************************** Submenu routine **************************************/
+
+/** Submenu for adding atoms to the molecule.
+ * \param *periode periodentafel
+ * \param *mol the molecule to add to
+ */
+void AddAtoms(periodentafel *periode, molecule *mol)
+{
+  atom *first, *second, *third, *fourth;
+  vector **atoms;
+  vector x,y,z,n;  // coordinates for absolute point in cell volume
+  double a,b,c;
+  char choice;  // menu choice char
+  bool valid;
+
+  cout << Verbose(0) << "===========ADD ATOM============================" << endl;
+  cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
+  cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
+  cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
+  cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
+  cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+  
+  switch (choice) {
+      case 'a': // absolute coordinates of atom
+        cout << Verbose(0) << "Enter absolute coordinates." << endl;
+        first = new atom;
+        first->x.AskPosition(mol->cell_size, false);
+        first->type = periode->AskElement();  // give type
+        mol->AddAtom(first);  // add to molecule
+        break;
+        
+      case 'b': // relative coordinates of atom wrt to reference point
+        first = new atom;
+        valid = true;
+        do {
+          if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
+          cout << Verbose(0) << "Enter reference coordinates." << endl;
+          x.AskPosition(mol->cell_size, true);
+          cout << Verbose(0) << "Enter relative coordinates." << endl;
+          first->x.AskPosition(mol->cell_size, false);
+          first->x.AddVector((const vector *)&x);
+          cout << Verbose(0) << "\n";
+        } while (!(valid = mol->CheckBounds((const vector *)&first->x)));
+        first->type = periode->AskElement();  // give type
+        mol->AddAtom(first);  // add to molecule
+        break;
+        
+      case 'c': // relative coordinates of atom wrt to already placed atom
+        first = new atom;
+        valid = true;
+        do {
+          if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
+          second = mol->AskAtom("Enter atom number: ");                
+          cout << Verbose(0) << "Enter relative coordinates." << endl;
+          first->x.AskPosition(mol->cell_size, false);
+          for (int i=0;i<3;i++) {
+            first->x.x[i] += second->x.x[i];
+          }
+        } while (!(valid = mol->CheckBounds((const vector *)&first->x)));
+        first->type = periode->AskElement();  // give type
+        mol->AddAtom(first);  // add to molecule
+        break;
+      
+      case 'd': // two atoms, two angles and a distance
+        first = new atom;
+        valid = true;
+        do {
+          if (!valid) {
+            cout << Verbose(0) << "Resulting coordinates out of cell - ";
+            first->x.Output((ofstream *)&cout);
+            cout << Verbose(0) << endl;
+          }
+          cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
+          second = mol->AskAtom("Enter central atom: ");
+          third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
+          fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
+          a = ask_value("Enter distance between central (first) and new atom: ");
+          b = ask_value("Enter angle between new, first and second atom (degrees): ");
+          b *= M_PI/180.;
+          bound(&b, 0., 2.*M_PI);
+          c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
+          c *= M_PI/180.;
+          bound(&c, -M_PI, M_PI);
+          cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
+/*
+          second->Output(1,1,(ofstream *)&cout);
+          third->Output(1,2,(ofstream *)&cout);
+          fourth->Output(1,3,(ofstream *)&cout);
+          n.MakeNormalVector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
+          x.CopyVector(&second->x);
+          x.SubtractVector(&third->x);
+          x.CopyVector(&fourth->x);
+          x.SubtractVector(&third->x);
+          
+          if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
+            cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
+            continue;
+          }
+          cout << Verbose(0) << "resulting relative coordinates: ";
+          z.Output((ofstream *)&cout);
+          cout << Verbose(0) << endl;
+          */
+          // calc axis vector
+          x.CopyVector(&second->x);
+          x.SubtractVector(&third->x);
+          x.Normalize();
+          cout << "x: ",
+          x.Output((ofstream *)&cout);
+          cout << endl;          
+          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
+          cout << "z: ",
+          z.Output((ofstream *)&cout);
+          cout << endl;          
+          y.MakeNormalVector(&x,&z);
+          cout << "y: ",
+          y.Output((ofstream *)&cout);
+          cout << endl;          
+         
+          // rotate vector around first angle
+          first->x.CopyVector(&x);
+          first->x.RotateVector(&z,b - M_PI);
+          cout << "Rotated vector: ",
+          first->x.Output((ofstream *)&cout);
+          cout << endl;          
+          // remove the projection onto the rotation plane of the second angle
+          n.CopyVector(&y);
+          n.Scale(first->x.Projection(&y));
+          cout << "N1: ",
+          n.Output((ofstream *)&cout);
+          cout << endl;          
+          first->x.SubtractVector(&n);
+          cout << "Subtracted vector: ",
+          first->x.Output((ofstream *)&cout);
+          cout << endl;          
+          n.CopyVector(&z);
+          n.Scale(first->x.Projection(&z));
+          cout << "N2: ",
+          n.Output((ofstream *)&cout);
+          cout << endl;          
+          first->x.SubtractVector(&n);
+          cout << "2nd subtracted vector: ",
+          first->x.Output((ofstream *)&cout);
+          cout << endl;          
+          
+          // rotate another vector around second angle
+          n.CopyVector(&y);
+          n.RotateVector(&x,c - M_PI);
+          cout << "2nd Rotated vector: ",
+          n.Output((ofstream *)&cout);
+          cout << endl;          
+          
+          // add the two linear independent vectors
+          first->x.AddVector(&n);
+          first->x.Normalize();         
+          first->x.Scale(a);
+          first->x.AddVector(&second->x);
+          
+          cout << Verbose(0) << "resulting coordinates: ";
+          first->x.Output((ofstream *)&cout);
+          cout << Verbose(0) << endl;
+        } while (!(valid = mol->CheckBounds((const vector *)&first->x)));
+        first->type = periode->AskElement();  // give type
+        mol->AddAtom(first);  // add to molecule
+        break;
+
+      case 'e': // least square distance position to a set of atoms
+        first = new atom;
+        atoms = new (vector*[128]);
+        valid = true;
+        for(int i=0;i<128;i++)
+          atoms[i] = NULL;
+        int i=0, j=0;
+        cout << Verbose(0) << "Now we need at least three molecules.\n";
+        do {
+          cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
+          cin >> j;
+          if (j != -1) {
+            second = mol->FindAtom(j);
+            atoms[i++] = &(second->x);
+          }
+        } while ((j != -1) && (i<128));
+        if (i >= 2) {
+          first->x.LSQdistance(atoms, i);             
+
+          first->x.Output((ofstream *)&cout);
+          first->type = periode->AskElement();  // give type
+          mol->AddAtom(first);  // add to molecule
+        } else {
+          delete first;
+          cout << Verbose(0) << "Please enter at least two vectors!\n";
+        }
+        break;
+  };
+};
+
+/** Submenu for centering the atoms in the molecule.
+ * \param *mol the molecule with all the atoms
+ */
+void CenterAtoms(molecule *mol)
+{
+  vector x;
+  char choice;  // menu choice char
+  double min[3];
+  int j;
+  
+  cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - on origin" << endl;
+  cout << Verbose(0) << " b - on center of gravity" << endl;
+  cout << Verbose(0) << " c - within box with additional boundary" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+  
+  switch (choice) {
+    default:
+      cout << Verbose(0) << "Not a valid choice." << endl;
+      break;
+    case 'a':
+      cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
+      mol->CenterOrigin((ofstream *)&cout, &x);
+      break;
+    case 'b':
+      cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
+      mol->CenterGravity((ofstream *)&cout, &x);
+      break;
+    case 'c':
+      cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
+      for (int i=0;i<3;i++) {
+        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
+        cin >> min[i];
+      }
+      mol->CenterEdge((ofstream *)&cout, &x);  // make every coordinate positive
+      mol->SetBoxDimension(&x);  // update Box of atoms by boundary
+      // translate each coordinate by boundary
+      j=-1;
+      for (int i=0;i<NDIM;i++) {
+        j += i+1;
+        x.x[i] = min[i];
+        mol->cell_size[j] += x.x[i]*2.;
+      }
+      mol->Translate(&x);
+      break;
+  }
+};
+
+/** Submenu for aligning the atoms in the molecule.
+ * \param *periode periodentafel
+ * \param *mol the molecule with all the atoms
+ */
+void AlignAtoms(periodentafel *periode, molecule *mol)
+{
+  atom *first, *second, *third;
+  vector x,n;
+  char choice;  // menu choice char
+
+  cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
+  cout << Verbose(0) << " b - state alignment vector" << endl;
+  cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
+  cout << Verbose(0) << " d - align automatically by least square fit" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+  
+  switch (choice) {
+    default:
+    case 'a': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+      third = mol->AskAtom("Enter third atom: ");
+
+      n.MakeNormalVector((const vector *)&first->x,(const vector *)&second->x,(const vector *)&third->x);
+      break;
+    case 'b': // normal vector of mirror plane
+      cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
+      n.AskPosition(mol->cell_size,0);
+      n.Normalize();
+      break;
+    case 'c': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+
+      n.CopyVector((const vector *)&first->x);  
+      n.SubtractVector((const vector *)&second->x);  
+      n.Normalize();
+      break;        
+    case 'd':
+      char shorthand[4];
+      vector a;
+      struct lsq_params param;
+      do {
+        fprintf(stdout, "Enter the element of atoms to be chosen: ");
+        fscanf(stdin, "%3s", shorthand);
+      } while ((param.type = periode->FindElement(shorthand)) == NULL);
+      cout << Verbose(0) << "Element is " << param.type->name << endl;
+      mol->GetAlignVector(&param);
+      for (int i=0;i<3;i++) {
+        x.x[i] = gsl_vector_get(param.x,i);
+        n.x[i] = gsl_vector_get(param.x,i+3);
+      } 
+      gsl_vector_free(param.x);
+      cout << Verbose(0) << "Offset vector: ";
+      x.Output((ofstream *)&cout);
+      cout << Verbose(0) << endl;
+      n.Normalize();
+      break;       
+  };
+  cout << Verbose(0) << "Alignment vector: ";
+  n.Output((ofstream *)&cout);
+  cout << Verbose(0) << endl;
+  mol->Align(&n);
+};
+
+/** Submenu for mirroring the atoms in the molecule.
+ * \param *mol the molecule with all the atoms
+ */
+void MirrorAtoms(molecule *mol)
+{
+  atom *first, *second, *third;
+  vector n;
+  char choice;  // menu choice char
+  
+  cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
+  cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
+  cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+  
+  switch (choice) {
+    default:
+    case 'a': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+      third = mol->AskAtom("Enter third atom: ");
+
+      n.MakeNormalVector((const vector *)&first->x,(const vector *)&second->x,(const vector *)&third->x);
+      break;
+    case 'b': // normal vector of mirror plane
+      cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
+      n.AskPosition(mol->cell_size,0);
+      n.Normalize();
+      break;
+    case 'c': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+
+      n.CopyVector((const vector *)&first->x);  
+      n.SubtractVector((const vector *)&second->x);  
+      n.Normalize();
+      break;          
+  };
+  cout << Verbose(0) << "Normal vector: ";
+  n.Output((ofstream *)&cout);
+  cout << Verbose(0) << endl;
+  mol->Mirror((const vector *)&n);
+};
+
+/** Submenu for removing the atoms from the molecule.
+ * \param *mol the molecule with all the atoms
+ */
+void RemoveAtoms(molecule *mol)
+{
+  atom *first, *second;
+  int axis;
+  double tmp1, tmp2;
+  char choice;  // menu choice char
+  
+  cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state atom for removal by number" << endl;
+  cout << Verbose(0) << " b - keep only in radius around atom" << endl;
+  cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+  
+  switch (choice) {
+    default:
+    case 'a':
+      if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
+        cout << Verbose(1) << "Atom removed." << endl;
+      else
+        cout << Verbose(1) << "Atom not found." << endl;
+      break;
+    case 'b':
+      second = mol->AskAtom("Enter number of atom as reference point: ");
+      cout << Verbose(0) << "Enter radius: ";
+      cin >> tmp1;
+      first = mol->start;
+      while(first->next != mol->end) {
+        first = first->next;
+        if (first->x.Distance((const vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
+          mol->RemoveAtom(first);
+      }
+      break;
+    case 'c':
+      cout << Verbose(0) << "Which axis is it: ";
+      cin >> axis;
+      cout << Verbose(0) << "Left inward boundary: ";
+      cin >> tmp1;
+      cout << Verbose(0) << "Right inward boundary: ";
+      cin >> tmp2;
+      first = mol->start;
+      while(first->next != mol->end) {
+        first = first->next;
+        if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...
+          mol->RemoveAtom(first);
+      }
+      break;          
+  };
+  //mol->Output((ofstream *)&cout);
+  choice = 'r';
+};
+
+/** Submenu for measuring out the atoms in the molecule.
+ * \param *periode periodentafel
+ * \param *mol the molecule with all the atoms
+ */
+void MeasureAtoms(periodentafel *periode, molecule *mol)
+{
+  atom *first, *second, *third;
+  vector x,y;
+  double min[256], tmp1, tmp2, tmp3;
+  int Z;
+  char choice;  // menu choice char
+  
+  cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
+  cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
+  cout << Verbose(0) << " c - calculate bond angle" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch(choice) {
+    default:
+      cout << Verbose(1) << "Not a valid choice." << endl;
+      break;
+    case 'a':
+      first = mol->AskAtom("Enter first atom: ");
+      for (int i=0;i<256;i++)
+        min[i] = 0.;
+        
+      second = mol->start;   
+      while ((second->next != mol->end)) {
+        second = second->next; // advance
+        Z = second->type->Z;
+        tmp1 = 0.;
+        if (first != second) {
+          x.CopyVector((const vector *)&first->x);
+          x.SubtractVector((const vector *)&second->x);
+          tmp1 = x.Norm();
+        }
+        if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
+        //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;         
+      }
+      for (int i=0;i<256;i++)
+        if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
+      break;
+      
+    case 'b':
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+      for (int i=0;i<NDIM;i++)
+        min[i] = 0.;
+      x.CopyVector((const vector *)&first->x);
+      x.SubtractVector((const vector *)&second->x);
+      tmp1 = x.Norm();
+      cout << Verbose(1) << "Distance vector is ";
+      x.Output((ofstream *)&cout);
+      cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
+      break;
+
+    case 'c':
+      cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter central atom: ");
+      third  = mol->AskAtom("Enter last atom: ");
+      tmp1 = tmp2 = tmp3 = 0.;
+      x.CopyVector((const vector *)&first->x);
+      x.SubtractVector((const vector *)&second->x);
+      y.CopyVector((const vector *)&third->x);
+      y.SubtractVector((const vector *)&second->x);
+      cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
+      cout << Verbose(0) << (acos(x.ScalarProduct((const vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;         
+      break;
+  }
+};
+
+/** Submenu for measuring out the atoms in the molecule.
+ * \param *mol the molecule with all the atoms
+ * \param *configuration configuration structure for the to be written config files of all fragments
+ */
+void FragmentAtoms(molecule *mol, config *configuration)
+{
+  enum BondOrderScheme Scheme = NoScheme;
+  enum CutCyclicBond CutCyclic;
+  char schema;
+  int Order1, Order2;
+  clock_t start, end;
+  
+  cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
+  cout << Verbose(0) << "What's the desired bond order: ";
+  cin >> Order1;
+  cout << Verbose(0) << "What's the desired bond scheme [(B)ottomUp/(T)opDown/(A)NOVA/(C)ombined]: ";
+  cin >> schema;
+  CutCyclic = SaturateBond;
+  switch (schema) {
+    default:
+        Scheme = NoScheme;
+      break;
+    case 'B':
+        CutCyclic = KeepBond;
+    case 'b':
+        Scheme = BottomUp;
+        Order2 = 0;
+      break;
+    case 'T':
+        CutCyclic = KeepBond;
+    case 't':
+        Scheme = TopDown;
+        Order2 = Order1;
+        Order1 = 0;
+      break;
+    case 'A':
+        CutCyclic = KeepBond;
+    case 'a':
+        Scheme = ANOVA;
+        Order2 = 0;
+      break;
+    case 'C':
+        CutCyclic = KeepBond;
+    case 'c':
+        cout << Verbose(0) << "Combined first selects subgraphs by BottomUp of bond order " << Order1 << " then TopDown to fragment these." << endl;
+        cout << Verbose(0) << "What's the desired bond order for TopDown: ";
+        cin >> Order2;
+        Scheme = Combined;
+      break;
+  };
+  if (mol->first->next != mol->last) {  // there are bonds
+    start = clock();
+    mol->FragmentMolecule((ofstream *)&cout, Order1, Order2, Scheme, configuration, CutCyclic);
+    end = clock();
+    cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+  } else 
+    cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
+};
+
+/********************************************** Test routine **************************************/
+
+/** Is called always as option 'T' in the menu.
+ */
+void testroutine(molecule *mol)
+{
+  // the current test routine checks the functionality of the KeySet&Graph concept:
+  // We want to have a multiindex (the KeySet) describing a unique subgraph
+  atom *Walker = mol->start;
+  int i, comp, counter=0;
+  
+  // generate some KeySets
+  cout << "Generating KeySets." << endl;
+  KeySet TestSets[mol->AtomCount+1];
+  i=1;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    for (int j=0;j<i;j++) {
+      TestSets[j].insert(Walker->nr);
+    }
+    i++;
+  }
+  cout << "Testing insertion of already present item in KeySets." << endl;
+  KeySetTestPair test;
+  test = TestSets[mol->AtomCount-1].insert(Walker->nr); 
+  if (test.second) {
+    cout << Verbose(1) << "Insertion worked?!" << endl;
+  } else {
+    cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
+  }
+  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
+  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
+
+  // constructing Graph structure  
+  cout << "Generating Subgraph class." << endl;
+  Graph Subgraphs;
+
+  // insert KeySets into Subgraphs
+  cout << "Inserting KeySets into Subgraph class." << endl;
+  for (int j=0;j<mol->AtomCount;j++) {
+    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
+  }
+  cout << "Testing insertion of already present item in Subgraph." << endl;
+  GraphTestPair test2;
+  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.))); 
+  if (test2.second) {
+    cout << Verbose(1) << "Insertion worked?!" << endl;
+  } else {
+    cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
+  }
+  
+  // show graphs
+  cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
+  Graph::iterator A = Subgraphs.begin();
+  while (A !=  Subgraphs.end()) {
+    cout << (*A).second.first << ": ";
+    KeySet::iterator key = (*A).first.begin();
+    comp = -1;
+    while (key != (*A).first.end()) {
+      if ((*key) > comp)
+        cout << (*key) << " ";
+      else 
+        cout << (*key) << "! ";
+      comp = (*key);
+      key++;
+    }
+    cout << endl;
+    A++;
+  }
+};
+
+/********************************************** Main routine **************************************/
+ 
+int main(int argc, char **argv)
+{
+  periodentafel *periode = new periodentafel; // and a period table of all elements
+  molecule *mol = new molecule(periode);    // first we need an empty molecule
+  double tmp1;
+  double bond, min_bond;
+  atom *first, *second;
+  element *finder;
+  char choice;  // menu choice char
+  vector x,y,z,n;  // coordinates for absolute point in cell volume
+	double *factor;	// unit factor if desired
+  bool valid; // flag if input was valid or not
+  config configuration;
+  ifstream test;
+  ofstream output;
+  string line;
+  char filename[255];
+  int flag = 1;
+  int Z;
+  int j, axis, count, faktor;
+  int MinimumRingSize = -1;
+  enum BondOrderScheme Scheme = NoScheme;
+  enum CutCyclicBond CutCyclic;
+  MoleculeLeafClass *Subgraphs = NULL;
+  clock_t start,end;
+  element **Elements;
+  vector **Vectors;
+
+  // Greetings
+  cout << Verbose(0) << "Molecuilder v1" << endl << endl;
+  cout << Verbose(0) << "Build your own molecule position set." << endl;
+
+  // load element list
+  periode->LoadPeriodentafel();
+  cout << Verbose(1) << "Element list loaded." << endl;  
+  
+  // parse config file if given
+  if (argc >= 2) { // config file specified as option
+  	cout << Verbose(0) << "Config file given." << endl;
+    test.open(argv[1], ios::in);
+    if (test == NULL) {
+      cout << Verbose(1) << "Specified config file not found." << endl;
+      //return (1);
+    } else {
+      cout << Verbose(1) << "Specified config file found, parsing ...";
+      switch (configuration.TestSyntax(&test, periode, mol)) {
+        case 1:
+          cout << "new syntax." << endl;
+          configuration.Load(&test, periode, mol);
+          break;
+        case 0:
+          cout << "old syntax." << endl;
+          configuration.LoadOld(&test, periode, mol);
+          break;
+        default:
+          cout << "unknown syntax or empty configuration file." << endl;
+          //return(1);
+      }
+      test.close();
+    }
+  }
+  // command line stuff
+  if (argc > 2) {
+    cout << Verbose(0) << "Recognized command line argument: " << argv[2][1] << ".\n";
+    switch(argv[2][1]) {
+      case 't':
+        cout << Verbose(1) << "Translating all ions to new origin." << endl;
+        for (int i=0;i<3;i++)
+          x.x[i] = atof(argv[3+i]);
+        mol->Translate((const vector *)&x);
+        break;
+      case 'a':
+        cout << Verbose(1) << "Adding new atom." << endl;
+        first = new atom;
+        for (int i=0;i<3;i++)
+          first->x.x[i] = atof(argv[4+i]);
+        finder = periode->start;
+        while (finder != periode->end) {
+          finder = finder->next;
+          if (strncmp(finder->symbol,argv[3],3) == 0) {
+            first->type = finder;
+            break;
+          }
+        }
+        mol->AddAtom(first);  // add to molecule
+        break;
+      case 's':
+        j = -1;
+        cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
+        factor = (double *) Malloc(sizeof(double)*NDIM, "main: *factor");
+        factor[0] = atof(argv[3]);
+        factor[1] = (argc >=4) ? atof(argv[4]) : atof(argv[3]); 
+        factor[2] = (argc >=5) ? atof(argv[5]) : atof(argv[3]); 
+        mol->Scale(&factor);
+        for (int i=0;i<3;i++) {
+          j += i+1;
+          x.x[i] = atof(argv[3+i]);
+          mol->cell_size[j]*=factor[i];
+        }
+        Free((void **)&factor, "main: *factor");
+        break;
+      case 'c':
+        j = -1;
+        cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
+        // make every coordinate positive
+        mol->CenterEdge((ofstream *)&cout, &x);
+        // update Box of atoms by boundary
+        mol->SetBoxDimension(&x);
+        // translate each coordinate by boundary
+        j=-1;
+        for (int i=0;i<3;i++) {
+          j += i+1;
+          x.x[i] = atof(argv[3+i]);
+          mol->cell_size[j] += x.x[i]*2.;
+        }
+        mol->Translate((const vector *)&x);
+        break;
+      case 'r':
+        cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
+        break;
+      case 'p':
+        cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
+        if (!mol->AddXYZFile(argv[3]))
+          cout << Verbose(2) << "File not found." << endl;
+        else
+          cout << Verbose(2) << "File found and parsed." << endl;
+        break;
+      case 'f':
+        int i,j;
+        flag = 0;
+        if (argc > 4) {
+          cout << Verbose(0) << "Creating connection matrix..." << endl;
+          start = clock();
+          mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[3]));
+          cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
+          cout << Verbose(0) << "What's the desired bond scheme [(B)ottomUp/(T)opDown/(A)NOVA/(C)ombined]: ";
+          CutCyclic = SaturateBond;
+          switch (argv[5][0]) {
+            default:
+                Scheme = NoScheme;
+              break;
+            case 'B':
+              CutCyclic = KeepBond;
+            case 'b':
+                Scheme = BottomUp;
+                i = atoi(argv[4]);
+                j = 0;
+              break;
+            case 'T':
+              CutCyclic = KeepBond;
+            case 't':
+                Scheme = TopDown;
+                i = 0;
+                j = atoi(argv[4]);
+              break;
+            case 'A':
+              CutCyclic = KeepBond;
+            case 'a':
+                Scheme = ANOVA;
+                i = atoi(argv[4]);
+                j = 0;
+              break;
+            case 'C':
+              CutCyclic = KeepBond;
+            case 'c':
+                if (argc > 5) {
+                  i = atoi(argv[4]);
+                  j = atoi(argv[6]);
+                  Scheme = Combined;
+                }
+                else
+                  cerr << "Missing second bond order for TopDown fragmentation in combined approach." << endl;
+              break;
+          };
+          if (mol->first->next != mol->last) {
+            mol->FragmentMolecule((ofstream *)&cout, i, j, Scheme, &configuration, CutCyclic);
+          }
+          end = clock();
+          cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+        }
+        break;                
+      default:
+        break;
+    }
+    if (flag) {
+      mol->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
+      mol->CalculateOrbitals(configuration);
+      configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
+      if (argc >= 1)
+        output.open(argv[1], ios::trunc);
+      if (!output)
+        output.open("main_pcp_linux", ios::trunc);
+      if (configuration.Save(&output, periode, mol))
+        cout << Verbose(1) << "Saving of config file successful." << endl;
+      else
+        cout << Verbose(1) << "Saving of config file failed." << endl;
+      output.close();
+      output.clear();
+      // and save to xyz file
+      if (argc >= 1) {
+        strcpy(filename, argv[1]);
+        strcat(filename, ".xyz");
+        output.open(filename, ios::trunc);
+      }
+      if (output == NULL) {
+        strcpy(filename,"main_pcp_linux");
+        strcat(filename, ".xyz");
+        output.open(filename, ios::trunc);
+      } 
+      if (mol->OutputXYZ(&output))
+        cout << Verbose(0) << "Saving of XYZ file successful." << endl;
+      else
+        cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+      output.close();
+      output.clear();
+    }
+    delete(mol);
+    delete(periode);
+    return (0);
+  }
+
+  
+  // General stuff
+  if (mol->cell_size[0] == 0.) {
+    cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
+    for (int i=0;i<6;i++) {
+      cout << Verbose(1) << "Cell size" << i << ": ";
+      cin >> mol->cell_size[i];
+    }
+  }
+
+  // now the main construction loop
+  cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
+  do {    
+    cout << Verbose(0) << endl << endl;
+    cout << Verbose(0) << "============Element list=======================" << endl;
+    mol->Checkout((ofstream *)&cout);
+    cout << Verbose(0) << "============Atom list==========================" << endl;
+    mol->Output((ofstream *)&cout);
+    cout << Verbose(0) << "============Menu===============================" << endl;
+    cout << Verbose(0) << "a - add an atom" << endl;
+    cout << Verbose(0) << "r - remove an atom" << endl;
+    cout << Verbose(0) << "b - scale a bond between atoms" << endl;
+    cout << Verbose(0) << "u - change an atoms element" << endl;
+    cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
+    cout << Verbose(0) << "-----------------------------------------------" << endl;
+    cout << Verbose(0) << "p - Parse xyz file" << endl;
+    cout << Verbose(0) << "e - edit the current configuration" << endl;
+    cout << Verbose(0) << "o - create connection matrix" << endl;
+    cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
+    cout << Verbose(0) << "-----------------------------------------------" << endl;
+    cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
+    cout << Verbose(0) << "i - realign molecule" << endl; 
+    cout << Verbose(0) << "m - mirror all molecules" << endl; 
+    cout << Verbose(0) << "t - translate molecule by vector" << endl;
+    cout << Verbose(0) << "c - scale by unit transformation" << endl;
+    cout << Verbose(0) << "g - center atoms in box" << endl;
+    cout << Verbose(0) << "-----------------------------------------------" << endl;
+    cout << Verbose(0) << "s - save current setup to config file" << endl;
+    cout << Verbose(0) << "T - call the current test routine" << endl;
+    cout << Verbose(0) << "q - quit" << endl;
+    cout << Verbose(0) << "===============================================" << endl;
+    cout << Verbose(0) << "Input: ";
+    cin >> choice;
+    
+    switch (choice) {
+      default:
+      case 'q': // quit 
+        break;
+        
+      case 'a': // add atom
+        AddAtoms(periode, mol);
+        choice = 'a'; 
+        break;
+      
+      case 'd': // duplicate the periodic cell along a given axis, given times
+        cout << Verbose(0) << "State the axis [(+-)123]: ";
+        cin >> axis;
+        cout << Verbose(0) << "State the factor: ";
+        cin >> faktor;
+        
+        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 = (element **) Malloc(sizeof(element *)*count, "main: duplicateCell - **Elements");
+          Vectors = (vector **) Malloc(sizeof(vector *)*count, "main: duplicateCell - **Vectors");
+          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=0;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, ...)
+            }
+          }
+          if (mol->first->next != mol->last) // if connect matrix is present already, redo it
+            mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance);
+          // free memory
+          Free((void **)&Elements, "main: duplicateCell - **Elements");
+          Free((void **)&Vectors, "main: duplicateCell - **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;
+      
+      case 'g': // center the atoms
+        CenterAtoms(mol);
+        break;
+        
+      case 'b': // scale a bond
+        cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
+        first = mol->AskAtom("Enter first (fixed) atom: ");
+        second = mol->AskAtom("Enter second (shifting) atom: ");
+        min_bond = 0.;
+        for (int i=0;i<3;i++)
+          min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
+        min_bond = sqrt(min_bond);
+        cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
+        cout << Verbose(0) << "Enter new bond length [a.u.]: ";
+        cin >> bond;
+        for (int i=0;i<3;i++) {
+          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
+        }
+        //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
+        //second->Output(second->type->No, 1, (ofstream *)&cout);       
+        break;
+
+      case 'i': // align all atoms 
+        AlignAtoms(periode, mol);
+        break;
+
+      case 'm': // mirror atoms along a given axis
+        MirrorAtoms(mol);
+        break;
+       
+      case 't': // translate all atoms
+	     cout << Verbose(0) << "Enter translation vector." << endl;       
+	     x.AskPosition(mol->cell_size,0);
+       mol->Translate((const vector *)&x);
+	     break;
+  
+      case 'e': // edit each field of the configuration
+       configuration.Edit(mol);
+       break;
+	
+      case 'c': // unit scaling of the metric 
+    	 cout << Verbose(0) << "Enter three factors: ";
+       factor = (double *) Malloc(sizeof(double)*NDIM, "main: *factor");
+       cin >> factor[0];
+       cin >> factor[1];
+       cin >> factor[2];
+       valid = true;
+       mol->Scale(&factor);
+	     Free((void **)&factor, "main: *factor");
+       break;
+			 
+      case 'r': // remove atom
+        RemoveAtoms(mol);        
+        break;
+        
+      case 'l': // measure distances or angles
+        MeasureAtoms(periode, mol);
+        break;
+
+      case 'p': // parse and XYZ file
+        cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
+        do {
+          cout << Verbose(0) << "Enter file name: ";
+          cin >> filename;
+        } while (!mol->AddXYZFile(filename));
+        break;
+
+      case 'o': // create the connection matrix
+        cout << Verbose(0) << "What's the maximum bond distance: ";
+        cin >> tmp1;
+        start = clock();
+        mol->CreateAdjacencyList((ofstream *)&cout, tmp1);
+        //mol->CreateListOfBondsPerAtom((ofstream *)&cout);
+        Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
+        while (Subgraphs->next != NULL) {
+          Subgraphs = Subgraphs->next;
+          delete(Subgraphs->previous);
+        }
+        delete(Subgraphs);    // we don't need the list here, so free everything
+        Subgraphs = NULL;
+        end = clock();
+        cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+        break;
+        
+      case 'f':
+        FragmentAtoms(mol, &configuration);
+      	break;
+        
+      case 'u': // change an atom's element
+        first = NULL;
+        do {
+          cout << Verbose(0) << "Change the element of which atom: ";
+          cin >> Z;
+        } while ((first = mol->FindAtom(Z)) == NULL);
+        cout << Verbose(0) << "New element by atomic number Z: ";        
+        cin >> Z;
+        first->type = periode->FindElement(Z);
+        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;    
+        break;
+        
+      case 'T':
+        testroutine(mol);
+        break;
+      
+      case 's': // save to config file
+        // get correct valence orbitals
+        mol->CalculateOrbitals(configuration);
+        configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
+        if (argc >= 1)
+          output.open(argv[1], ios::trunc);
+        if (output == NULL)
+          output.open("main_pcp_linux", ios::trunc);
+        if (configuration.Save(&output, periode, mol))
+          cout << Verbose(0) << "Saving of config file successful." << endl;
+        else
+          cout << Verbose(0) << "Saving of config file failed." << endl;
+        output.close();
+        output.clear();
+        // and save to xyz file
+        if (argc >= 1) {
+          strcpy(filename, argv[1]);
+          strcat(filename, ".xyz");
+          output.open(filename, ios::trunc);
+        }
+        if (output == NULL) {
+          strcpy(filename,"main_pcp_linux");
+          strcat(filename, ".xyz");
+          output.open(filename, ios::trunc);
+        } 
+        if (mol->OutputXYZ(&output))
+          cout << Verbose(0) << "Saving of XYZ file successful." << endl;
+        else
+          cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+        output.close();
+        output.clear();
+        break;
+    };
+  } while (choice != 'q');
+  
+  // save element data base
+  if (periode->StorePeriodentafel())
+    cout << Verbose(0) << "Saving of elements.db successful." << endl;
+  else
+    cout << Verbose(0) << "Saving of elements.db failed." << endl;
+
+  // Free all
+  if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
+    while (Subgraphs->next != NULL) {
+      Subgraphs = Subgraphs->next;
+      delete(Subgraphs->previous);
+    }
+    delete(Subgraphs);
+  }
+  delete(mol);
+  delete(periode);
+  return (0);
+}
+
+/********************************************** E N D **************************************************/
