source: src/builder.cpp@ b60a29

/** \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 <cstring>

#include "analysis_bonds.hpp"
#include "analysis_correlation.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "bondgraph.hpp"
#include "boundary.hpp"
#include "config.hpp"
#include "element.hpp"
#include "ellipsoid.hpp"
#include "helpers.hpp"
#include "leastsquaremin.hpp"
#include "linkedcell.hpp"
#include "log.hpp"
#include "memoryusageobserverunittest.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "version.h"
#include "World.hpp"

/********************************************* Subsubmenu routine ************************************/

/** Submenu for adding atoms to the molecule.
 * \param *periode periodentafel
 * \param *molecule molecules with atoms
 */
static 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) {
    default:
      DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl);
      break;
      case 'a': // absolute coordinates of atom
        cout << Verbose(0) << "Enter absolute coordinates." << endl;
        first = new atom;
        first->x.AskPosition(World::get()->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) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
          cout << Verbose(0) << "Enter reference coordinates." << endl;
          x.AskPosition(World::get()->cell_size, true);
          cout << Verbose(0) << "Enter relative coordinates." << endl;
          first->x.AskPosition(World::get()->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) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
          second = mol->AskAtom("Enter atom number: ");
          Log() << Verbose(0) << "Enter relative coordinates." << endl;
          first->x.AskPosition(World::get()->cell_size, false);
          for (int i=NDIM;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) {
            DoeLog(2) && (eLog()<< Verbose(2) << "Resulting coordinates out of cell - " << first->x << 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)) {
         coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
            continue;
          }
          Log() << Verbose(0) << "resulting relative coordinates: ";
          z.Output();
          Log() << Verbose(0) << endl;
          */
          // calc axis vector
          x.CopyVector(&second->x);
          x.SubtractVector(&third->x);
          x.Normalize();
          Log() << Verbose(0) << "x: ",
          x.Output();
          Log() << Verbose(0) << endl;
          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
          Log() << Verbose(0) << "z: ",
          z.Output();
          Log() << Verbose(0) << endl;
          y.MakeNormalVector(&x,&z);
          Log() << Verbose(0) << "y: ",
          y.Output();
          Log() << Verbose(0) << endl;

          // rotate vector around first angle
          first->x.CopyVector(&x);
          first->x.RotateVector(&z,b - M_PI);
          Log() << Verbose(0) << "Rotated vector: ",
          first->x.Output();
          Log() << Verbose(0) << endl;
          // remove the projection onto the rotation plane of the second angle
          n.CopyVector(&y);
          n.Scale(first->x.ScalarProduct(&y));
          Log() << Verbose(0) << "N1: ",
          n.Output();
          Log() << Verbose(0) << endl;
          first->x.SubtractVector(&n);
          Log() << Verbose(0) << "Subtracted vector: ",
          first->x.Output();
          Log() << Verbose(0) << endl;
          n.CopyVector(&z);
          n.Scale(first->x.ScalarProduct(&z));
          Log() << Verbose(0) << "N2: ",
          n.Output();
          Log() << Verbose(0) << endl;
          first->x.SubtractVector(&n);
          Log() << Verbose(0) << "2nd subtracted vector: ",
          first->x.Output();
          Log() << Verbose(0) << endl;

          // rotate another vector around second angle
          n.CopyVector(&y);
          n.RotateVector(&x,c - M_PI);
          Log() << Verbose(0) << "2nd Rotated vector: ",
          n.Output();
          Log() << Verbose(0) << endl;

          // add the two linear independent vectors
          first->x.AddVector(&n);
          first->x.Normalize();
          first->x.Scale(a);
          first->x.AddVector(&second->x);

          Log() << Verbose(0) << "resulting coordinates: ";
          first->x.Output();
          Log() << 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=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((const Vector **)atoms, i);

          first->x.Output();
          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 molecule with all the atoms
 */
static void CenterAtoms(molecule *mol)
{
  Vector x, y, helper;
  char choice;  // menu choice char

  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) << " d - within given simulation box" << 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();
      break;
    case 'b':
      cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
      mol->CenterPeriodic();
      break;
    case 'c':
      cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
      for (int i=0;i<NDIM;i++) {
        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
        cin >> y.x[i];
      }
      mol->CenterEdge(&x);  // make every coordinate positive
      mol->Center.AddVector(&y); // translate by boundary
      helper.CopyVector(&y);
      helper.Scale(2.);
      helper.AddVector(&x);
      mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
      break;
    case 'd':
      cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
      for (int i=0;i<NDIM;i++) {
        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
        cin >> x.x[i];
      }
      // update Box of atoms by boundary
      mol->SetBoxDimension(&x);
      // center
      mol->CenterInBox();
      break;
  }
};

/** Submenu for aligning the atoms in the molecule.
 * \param *periode periodentafel
 * \param *mol molecule with all the atoms
 */
static 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(World::get()->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=NDIM;i--;) {
        x.x[i] = gsl_vector_get(param.x,i);
        n.x[i] = gsl_vector_get(param.x,i+NDIM);
      }
      gsl_vector_free(param.x);
      cout << Verbose(0) << "Offset vector: ";
      x.Output();
      Log() << Verbose(0) << endl;
      n.Normalize();
      break;
  };
  Log() << Verbose(0) << "Alignment vector: ";
  n.Output();
  Log() << Verbose(0) << endl;
  mol->Align(&n);
};

/** Submenu for mirroring the atoms in the molecule.
 * \param *mol molecule with all the atoms
 */
static void MirrorAtoms(molecule *mol)
{
  atom *first, *second, *third;
  Vector n;
  char choice;  // menu choice char

  Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
  Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
  Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl;
  Log() << Verbose(0) << " c - state two atoms in normal direction" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  Log() << 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
      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
      n.AskPosition(World::get()->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;
  };
  Log() << Verbose(0) << "Normal vector: ";
  n.Output();
  Log() << Verbose(0) << endl;
  mol->Mirror((const Vector *)&n);
};

/** Submenu for removing the atoms from the molecule.
 * \param *mol molecule with all the atoms
 */
static void RemoveAtoms(molecule *mol)
{
  atom *first, *second;
  int axis;
  double tmp1, tmp2;
  char choice;  // menu choice char

  Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
  Log() << Verbose(0) << " a - state atom for removal by number" << endl;
  Log() << Verbose(0) << " b - keep only in radius around atom" << endl;
  Log() << Verbose(0) << " c - remove this with one axis greater value" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch (choice) {
    default:
    case 'a':
      if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
        Log() << Verbose(1) << "Atom removed." << endl;
      else
        Log() << Verbose(1) << "Atom not found." << endl;
      break;
    case 'b':
      second = mol->AskAtom("Enter number of atom as reference point: ");
      Log() << Verbose(0) << "Enter radius: ";
      cin >> tmp1;
      first = mol->start;
      second = first->next;
      while(second != mol->end) {
        first = second;
        second = first->next;
        if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
          mol->RemoveAtom(first);
      }
      break;
    case 'c':
      Log() << Verbose(0) << "Which axis is it: ";
      cin >> axis;
      Log() << Verbose(0) << "Lower boundary: ";
      cin >> tmp1;
      Log() << Verbose(0) << "Upper boundary: ";
      cin >> tmp2;
      first = mol->start;
      second = first->next;
      while(second != mol->end) {
        first = second;
        second = first->next;
        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
          //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
          mol->RemoveAtom(first);
        }
      }
      break;
  };
  //mol->Output();
  choice = 'r';
};

/** Submenu for measuring out the atoms in the molecule.
 * \param *periode periodentafel
 * \param *mol molecule with all the atoms
 */
static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
{
  atom *first, *second, *third;
  Vector x,y;
  double min[256], tmp1, tmp2, tmp3;
  int Z;
  char choice;  // menu choice char

  Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
  Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
  Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl;
  Log() << Verbose(0) << " c - calculate bond angle" << endl;
  Log() << Verbose(0) << " d - calculate principal axis of the system" << endl;
  Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
  Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;
  Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch(choice) {
    default:
      Log() << Verbose(1) << "Not a valid choice." << endl;
      break;
    case 'a':
      first = mol->AskAtom("Enter first atom: ");
      for (int i=MAX_ELEMENTS;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;
        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
      }
      for (int i=MAX_ELEMENTS;i--;)
        if (min[i] != 0.) Log() << 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=NDIM;i--;)
        min[i] = 0.;
      x.CopyVector((const Vector *)&first->x);
      x.SubtractVector((const Vector *)&second->x);
      tmp1 = x.Norm();
      Log() << Verbose(1) << "Distance vector is ";
      x.Output();
      Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
      break;

    case 'c':
      Log() << 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);
      Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
      Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
      break;
    case 'd':
      Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;
      Log() << Verbose(0) << "Shall we rotate? [0/1]: ";
      cin >> Z;
      if ((Z >=0) && (Z <=1))
        mol->PrincipalAxisSystem((bool)Z);
      else
        mol->PrincipalAxisSystem(false);
      break;
    case 'e':
      {
        Log() << Verbose(0) << "Evaluating volume of the convex envelope.";
        class Tesselation *TesselStruct = NULL;
        const LinkedCell *LCList = NULL;
        LCList = new LinkedCell(mol, 10.);
        FindConvexBorder(mol, TesselStruct, LCList, NULL);
        double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
        Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;\
        delete(LCList);
        delete(TesselStruct);
      }
      break;
    case 'f':
      mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps);
      break;
    case 'g':
      {
        char filename[255];
        Log() << Verbose(0) << "Please enter filename: " << endl;
        cin >> filename;
        Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
        ofstream *output = new ofstream(filename, ios::trunc);
        if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
          Log() << Verbose(2) << "File could not be written." << endl;
        else
          Log() << Verbose(2) << "File stored." << endl;
        output->close();
        delete(output);
      }
      break;
  }
};

/** Submenu for measuring out the atoms in the molecule.
 * \param *mol molecule with all the atoms
 * \param *configuration configuration structure for the to be written config files of all fragments
 */
static void FragmentAtoms(molecule *mol, config *configuration)
{
  int Order1;
  clock_t start, end;

  Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
  Log() << Verbose(0) << "What's the desired bond order: ";
  cin >> Order1;
  if (mol->first->next != mol->last) {  // there are bonds
    start = clock();
    mol->FragmentMolecule(Order1, configuration);
    end = clock();
    Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
  } else
    Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
};

/********************************************** Submenu routine **************************************/

/** Submenu for manipulating atoms.
 * \param *periode periodentafel
 * \param *molecules list of molecules whose atoms are to be manipulated
 */
static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
{
  atom *first, *second, *third;
  molecule *mol = NULL;
  Vector x,y,z,n; // coordinates for absolute point in cell volume
  double *factor; // unit factor if desired
  double bond, minBond;
  char choice;  // menu choice char
  bool valid;

  Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl;
  Log() << Verbose(0) << "a - add an atom" << endl;
  Log() << Verbose(0) << "r - remove an atom" << endl;
  Log() << Verbose(0) << "b - scale a bond between atoms" << endl;
  Log() << Verbose(0) << "t - turn an atom round another bond" << endl;
  Log() << Verbose(0) << "u - change an atoms element" << endl;
  Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  if (molecules->NumberOfActiveMolecules() > 1)
    DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch (choice) {
    default:
      Log() << Verbose(0) << "Not a valid choice." << endl;
      break;

    case 'a': // add atom
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        AddAtoms(periode, mol);
      }
      break;

    case 'b': // scale a bond
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        Log() << Verbose(0) << "Scaling bond length between two atoms." << endl;
        first = mol->AskAtom("Enter first (fixed) atom: ");
        second = mol->AskAtom("Enter second (shifting) atom: ");
        minBond = 0.;
        for (int i=NDIM;i--;)
          minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
        minBond = sqrt(minBond);
        Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;
        Log() << Verbose(0) << "Enter new bond length [a.u.]: ";
        cin >> bond;
        for (int i=NDIM;i--;) {
          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
        }
        //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
        //second->Output(second->type->No, 1);
      }
      break;

    case 'c': // unit scaling of the metric
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
       Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
       Log() << Verbose(0) << "Enter three factors: ";
       factor = new double[NDIM];
       cin >> factor[0];
       cin >> factor[1];
       cin >> factor[2];
       valid = true;
       mol->Scale((const double ** const)&factor);
       delete[](factor);
      }
     break;

    case 'l': // measure distances or angles
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        MeasureAtoms(periode, mol, configuration);
      }
      break;

    case 'r': // remove atom
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        RemoveAtoms(mol);
      }
      break;

    case 't': // turn/rotate atom
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
          mol = *ListRunner;
          Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl;
          first = mol->AskAtom("Enter turning atom: ");
          second = mol->AskAtom("Enter central atom: ");
          third  = mol->AskAtom("Enter bond atom: ");
          cout << Verbose(0) << "Enter new angle in degrees: ";
          double tmp = 0.;
          cin >> tmp;
          // calculate old angle
          x.CopyVector((const Vector *)&first->x);
          x.SubtractVector((const Vector *)&second->x);
          y.CopyVector((const Vector *)&third->x);
          y.SubtractVector((const Vector *)&second->x);
          double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
          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) << alpha << " degrees" << endl;
          // rotate
          z.MakeNormalVector(&x,&y);
          x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
          x.AddVector(&second->x);
          first->x.CopyVector(&x);
          // check new angle
          x.CopyVector((const Vector *)&first->x);
          x.SubtractVector((const Vector *)&second->x);
          alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
          cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
          cout << Verbose(0) << alpha << " degrees" << endl;
        }
      break;

    case 'u': // change an atom's element
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        int Z;
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        first = NULL;
        do {
          Log() << Verbose(0) << "Change the element of which atom: ";
          cin >> Z;
        } while ((first = mol->FindAtom(Z)) == NULL);
        Log() << Verbose(0) << "New element by atomic number Z: ";
        cin >> Z;
        first->type = periode->FindElement(Z);
        Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
      }
      break;
  }
};

/** Submenu for manipulating molecules.
 * \param *periode periodentafel
 * \param *molecules list of molecule to manipulate
 */
static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
{
  atom *first = NULL;
  Vector x,y,z,n; // coordinates for absolute point in cell volume
  int j, axis, count, faktor;
  char choice;  // menu choice char
  molecule *mol = NULL;
  element **Elements;
  Vector **vectors;
  MoleculeLeafClass *Subgraphs = NULL;

  Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;
  Log() << Verbose(0) << "c - scale by unit transformation" << endl;
  Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
  Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
  Log() << Verbose(0) << "g - center atoms in box" << endl;
  Log() << Verbose(0) << "i - realign molecule" << endl;
  Log() << Verbose(0) << "m - mirror all molecules" << endl;
  Log() << Verbose(0) << "o - create connection matrix" << endl;
  Log() << Verbose(0) << "t - translate molecule by vector" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  if (molecules->NumberOfActiveMolecules() > 1)
    DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch (choice) {
    default:
      Log() << Verbose(0) << "Not a valid choice." << endl;
      break;

    case 'd': // duplicate the periodic cell along a given axis, given times
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        Log() << Verbose(0) << "State the axis [(+-)123]: ";
        cin >> axis;
        Log() << Verbose(0) << "State the factor: ";
        cin >> faktor;

        mol->CountAtoms(); // 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)
            DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
          x.Zero();
          y.Zero();
          y.x[abs(axis)-1] = World::get()->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, ...)
            }
          }
          if (mol->first->next != mol->last) // if connect matrix is present already, redo it
            mol->CreateAdjacencyList(mol->BondDistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
          // 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);
          }
          World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
        }
      }
      break;

    case 'f':
      FragmentAtoms(mol, configuration);
      break;

    case 'g': // center the atoms
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        CenterAtoms(mol);
      }
      break;

    case 'i': // align all atoms
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        AlignAtoms(periode, mol);
      }
      break;

    case 'm': // mirror atoms along a given axis
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        MirrorAtoms(mol);
      }
      break;

    case 'o': // create the connection matrix
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
          mol = *ListRunner;
          double bonddistance;
          clock_t start,end;
          Log() << Verbose(0) << "What's the maximum bond distance: ";
          cin >> bonddistance;
          start = clock();
          mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
          end = clock();
          Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
        }
      break;

    case 't': // translate all atoms
      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if ((*ListRunner)->ActiveFlag) {
        mol = *ListRunner;
        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
        Log() << Verbose(0) << "Enter translation vector." << endl;
        x.AskPosition(World::get()->cell_size,0);
        mol->Center.AddVector((const Vector *)&x);
     }
     break;
  }
  // 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);
  }
};


/** Submenu for creating new molecules.
 * \param *periode periodentafel
 * \param *molecules list of molecules to add to
 */
static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
{
  char choice;  // menu choice char
  Vector center;
  int nr, count;
  molecule *mol = NULL;

  Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl;
  Log() << Verbose(0) << "c - create new molecule" << endl;
  Log() << Verbose(0) << "l - load molecule from xyz file" << endl;
  Log() << Verbose(0) << "n - change molecule's name" << endl;
  Log() << Verbose(0) << "N - give molecules filename" << endl;
  Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl;
  Log() << Verbose(0) << "r - remove a molecule" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch (choice) {
    default:
      Log() << Verbose(0) << "Not a valid choice." << endl;
      break;
    case 'c':
      mol = new molecule(periode);
      molecules->insert(mol);
      break;

    case 'l': // load from XYZ file
      {
        char filename[MAXSTRINGSIZE];
        Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
        mol = new molecule(periode);
        do {
          Log() << Verbose(0) << "Enter file name: ";
          cin >> filename;
        } while (!mol->AddXYZFile(filename));
        mol->SetNameFromFilename(filename);
        // center at set box dimensions
        mol->CenterEdge(&center);
        double * const cell_size = World::get()->cell_size;
        cell_size[0] = center.x[0];
        cell_size[1] = 0;
        cell_size[2] = center.x[1];
        cell_size[3] = 0;
        cell_size[4] = 0;
        cell_size[5] = center.x[2];
        molecules->insert(mol);
      }
      break;

    case 'n':
      {
        char filename[MAXSTRINGSIZE];
        do {
          Log() << Verbose(0) << "Enter index of molecule: ";
          cin >> nr;
          mol = molecules->ReturnIndex(nr);
        } while (mol == NULL);
        Log() << Verbose(0) << "Enter name: ";
        cin >> filename;
        strcpy(mol->name, filename);
      }
      break;

    case 'N':
      {
        char filename[MAXSTRINGSIZE];
        do {
          Log() << Verbose(0) << "Enter index of molecule: ";
          cin >> nr;
          mol = molecules->ReturnIndex(nr);
        } while (mol == NULL);
        Log() << Verbose(0) << "Enter name: ";
        cin >> filename;
        mol->SetNameFromFilename(filename);
      }
      break;

    case 'p': // parse XYZ file
      {
        char filename[MAXSTRINGSIZE];
        mol = NULL;
        do {
          Log() << Verbose(0) << "Enter index of molecule: ";
          cin >> nr;
          mol = molecules->ReturnIndex(nr);
        } while (mol == NULL);
        Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
        do {
          Log() << Verbose(0) << "Enter file name: ";
          cin >> filename;
        } while (!mol->AddXYZFile(filename));
        mol->SetNameFromFilename(filename);
      }
      break;

    case 'r':
      Log() << Verbose(0) << "Enter index of molecule: ";
      cin >> nr;
      count = 1;
      for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
        if (nr == (*ListRunner)->IndexNr) {
          mol = *ListRunner;
          molecules->ListOfMolecules.erase(ListRunner);
          delete(mol);
          break;
        }
      break;
  }
};


/** Submenu for merging molecules.
 * \param *periode periodentafel
 * \param *molecules list of molecules to add to
 */
static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
{
  char choice;  // menu choice char

  Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl;
  Log() << Verbose(0) << "a - simple add of one molecule to another" << endl;
  Log() << Verbose(0) << "b - count the number of bonds of two elements" << endl;
  Log() << Verbose(0) << "B - count the number of bonds of three elements " << endl;
  Log() << Verbose(0) << "e - embedding merge of two molecules" << endl;
  Log() << Verbose(0) << "h - count the number of hydrogen bonds" << endl;
  Log() << Verbose(0) << "b - count the number of hydrogen bonds" << endl;
  Log() << Verbose(0) << "m - multi-merge of all molecules" << endl;
  Log() << Verbose(0) << "s - scatter merge of two molecules" << endl;
  Log() << Verbose(0) << "t - simple merge of two molecules" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch (choice) {
    default:
      Log() << Verbose(0) << "Not a valid choice." << endl;
      break;

    case 'a':
      {
        int src, dest;
        molecule *srcmol = NULL, *destmol = NULL;
        {
          do {
            Log() << Verbose(0) << "Enter index of destination molecule: ";
            cin >> dest;
            destmol = molecules->ReturnIndex(dest);
          } while ((destmol == NULL) && (dest != -1));
          do {
            Log() << Verbose(0) << "Enter index of source molecule to add from: ";
            cin >> src;
            srcmol = molecules->ReturnIndex(src);
          } while ((srcmol == NULL) && (src != -1));
          if ((src != -1) && (dest != -1))
            molecules->SimpleAdd(srcmol, destmol);
        }
      }
      break;

    case 'b':
      {
        const int nr = 2;
        char *names[nr] = {"first", "second"};
        int Z[nr];
        element *elements[nr];
        for (int i=0;i<nr;i++) {
          Z[i] = 0;
          do {
            cout << "Enter " << names[i] << " element: ";
            cin >> Z[i];
          } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
          elements[i] = periode->FindElement(Z[i]);
        }
        const int count = CountBondsOfTwo(molecules, elements[0], elements[1]);
        cout << endl << "There are " << count << " ";
        for (int i=0;i<nr;i++) {
          if (i==0)
            cout << elements[i]->symbol;
          else
            cout << "-" << elements[i]->symbol;
        }
        cout << " bonds." << endl;
      }
    break;

    case 'B':
      {
        const int nr = 3;
        char *names[nr] = {"first", "second", "third"};
        int Z[nr];
        element *elements[nr];
        for (int i=0;i<nr;i++) {
          Z[i] = 0;
          do {
            cout << "Enter " << names[i] << " element: ";
            cin >> Z[i];
          } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
          elements[i] = periode->FindElement(Z[i]);
        }
        const int count = CountBondsOfThree(molecules, elements[0], elements[1], elements[2]);
        cout << endl << "There are " << count << " ";
        for (int i=0;i<nr;i++) {
          if (i==0)
            cout << elements[i]->symbol;
          else
            cout << "-" << elements[i]->symbol;
        }
        cout << " bonds." << endl;
      }
    break;

    case 'e':
      {
        int src, dest;
        molecule *srcmol = NULL, *destmol = NULL;
        do {
          Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ";
          cin >> src;
          srcmol = molecules->ReturnIndex(src);
        } while ((srcmol == NULL) && (src != -1));
        do {
          Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ";
          cin >> dest;
          destmol = molecules->ReturnIndex(dest);
        } while ((destmol == NULL) && (dest != -1));
        if ((src != -1) && (dest != -1))
          molecules->EmbedMerge(destmol, srcmol);
      }
      break;

    case 'h':
      {
        int Z;
        cout << "Please enter interface element: ";
        cin >> Z;
        element * const InterfaceElement = periode->FindElement(Z);
        cout << endl << "There are " << CountHydrogenBridgeBonds(molecules, InterfaceElement) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << "." << endl;
      }
      break;

    case 'm':
      {
        int nr;
        molecule *mol = NULL;
        do {
          Log() << Verbose(0) << "Enter index of molecule to merge into: ";
          cin >> nr;
          mol = molecules->ReturnIndex(nr);
        } while ((mol == NULL) && (nr != -1));
        if (nr != -1) {
          int N = molecules->ListOfMolecules.size()-1;
          int *src = new int(N);
          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
            if ((*ListRunner)->IndexNr != nr)
              src[N++] = (*ListRunner)->IndexNr;        
          molecules->SimpleMultiMerge(mol, src, N);
          delete[](src);
        }
      }
      break;

    case 's':
      Log() << Verbose(0) << "Not implemented yet." << endl;
      break;

    case 't':
      {
        int src, dest;
        molecule *srcmol = NULL, *destmol = NULL;
        {
          do {
            Log() << Verbose(0) << "Enter index of destination molecule: ";
            cin >> dest;
            destmol = molecules->ReturnIndex(dest);
          } while ((destmol == NULL) && (dest != -1));
          do {
            Log() << Verbose(0) << "Enter index of source molecule to merge into: ";
            cin >> src;
            srcmol = molecules->ReturnIndex(src);
          } while ((srcmol == NULL) && (src != -1));
          if ((src != -1) && (dest != -1))
            molecules->SimpleMerge(srcmol, destmol);
        }
      }
      break;
  }
};


/********************************************** Test routine **************************************/

/** Is called always as option 'T' in the menu.
 * \param *molecules list of molecules
 */
static void testroutine(MoleculeListClass *molecules)
{
  // the current test routine checks the functionality of the KeySet&Graph concept:
  // We want to have a multiindex (the KeySet) describing a unique subgraph
  int i, comp, counter=0;

  // create a clone
  molecule *mol = NULL;
  if (molecules->ListOfMolecules.size() != 0) // clone
    mol = (molecules->ListOfMolecules.front())->CopyMolecule();
  else {
    DoeLog(0) && (eLog()<< Verbose(0) << "I don't have anything to test on ... ");
    performCriticalExit();
    return;
  }
  atom *Walker = mol->start;

  // generate some KeySets
  Log() << Verbose(0) << "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++;
  }
  Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
  KeySetTestPair test;
  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
  if (test.second) {
    Log() << Verbose(1) << "Insertion worked?!" << endl;
  } else {
    Log() << 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
  Log() << Verbose(0) << "Generating Subgraph class." << endl;
  Graph Subgraphs;

  // insert KeySets into Subgraphs
  Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl;
  for (int j=0;j<mol->AtomCount;j++) {
    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
  }
  Log() << Verbose(0) << "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) {
    Log() << Verbose(1) << "Insertion worked?!" << endl;
  } else {
    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
  }

  // show graphs
  Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl;
  Graph::iterator A = Subgraphs.begin();
  while (A !=  Subgraphs.end()) {
    Log() << Verbose(0) << (*A).second.first << ": ";
    KeySet::iterator key = (*A).first.begin();
    comp = -1;
    while (key != (*A).first.end()) {
      if ((*key) > comp)
        Log() << Verbose(0) << (*key) << " ";
      else
        Log() << Verbose(0) << (*key) << "! ";
      comp = (*key);
      key++;
    }
    Log() << Verbose(0) << endl;
    A++;
  }
  delete(mol);
};

/** Tries given filename or standard on saving the config file.
 * \param *ConfigFileName name of file
 * \param *configuration pointer to configuration structure with all the values
 * \param *periode pointer to periodentafel structure with all the elements
 * \param *molecules list of molecules structure with all the atoms and coordinates
 */
static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
{
  char filename[MAXSTRINGSIZE];
  ofstream output;
  molecule *mol = new molecule(periode);
  mol->SetNameFromFilename(ConfigFileName);

  if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
  }


  // first save as PDB data
  if (ConfigFileName != NULL)
    strcpy(filename, ConfigFileName);
  if (output == NULL)
    strcpy(filename,"main_pcp_linux");
  Log() << Verbose(0) << "Saving as pdb input ";
  if (configuration->SavePDB(filename, molecules))
    Log() << Verbose(0) << "done." << endl;
  else
    Log() << Verbose(0) << "failed." << endl;

  // then save as tremolo data file
  if (ConfigFileName != NULL)
    strcpy(filename, ConfigFileName);
  if (output == NULL)
    strcpy(filename,"main_pcp_linux");
  Log() << Verbose(0) << "Saving as tremolo data input ";
  if (configuration->SaveTREMOLO(filename, molecules))
    Log() << Verbose(0) << "done." << endl;
  else
    Log() << Verbose(0) << "failed." << endl;

  // translate each to its center and merge all molecules in MoleculeListClass into this molecule
  int N = molecules->ListOfMolecules.size();
  int *src = new int[N];
  N=0;
  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
    src[N++] = (*ListRunner)->IndexNr;
    (*ListRunner)->Translate(&(*ListRunner)->Center);
  }
  molecules->SimpleMultiAdd(mol, src, N);
  delete[](src);

  // ... and translate back
  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
    (*ListRunner)->Center.Scale(-1.);
    (*ListRunner)->Translate(&(*ListRunner)->Center);
    (*ListRunner)->Center.Scale(-1.);
  }

  Log() << Verbose(0) << "Storing configuration ... " << endl;
  // get correct valence orbitals
  mol->CalculateOrbitals(*configuration);
  configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
  if (ConfigFileName != NULL) { // test the file name
    strcpy(filename, ConfigFileName);
    output.open(filename, ios::trunc);
  } else if (strlen(configuration->configname) != 0) {
    strcpy(filename, configuration->configname);
    output.open(configuration->configname, ios::trunc);
    } else {
      strcpy(filename, DEFAULTCONFIG);
      output.open(DEFAULTCONFIG, ios::trunc);
    }
  output.close();
  output.clear();
  Log() << Verbose(0) << "Saving of config file ";
  if (configuration->Save(filename, periode, mol))
    Log() << Verbose(0) << "successful." << endl;
  else
    Log() << Verbose(0) << "failed." << endl;

  // and save to xyz file
  if (ConfigFileName != NULL) {
    strcpy(filename, ConfigFileName);
    strcat(filename, ".xyz");
    output.open(filename, ios::trunc);
  }
  if (output == NULL) {
    strcpy(filename,"main_pcp_linux");
    strcat(filename, ".xyz");
    output.open(filename, ios::trunc);
  }
  Log() << Verbose(0) << "Saving of XYZ file ";
  if (mol->MDSteps <= 1) {
    if (mol->OutputXYZ(&output))
      Log() << Verbose(0) << "successful." << endl;
    else
      Log() << Verbose(0) << "failed." << endl;
  } else {
    if (mol->OutputTrajectoriesXYZ(&output))
      Log() << Verbose(0) << "successful." << endl;
    else
      Log() << Verbose(0) << "failed." << endl;
  }
  output.close();
  output.clear();

  // and save as MPQC configuration
  if (ConfigFileName != NULL)
    strcpy(filename, ConfigFileName);
  if (output == NULL)
    strcpy(filename,"main_pcp_linux");
  Log() << Verbose(0) << "Saving as mpqc input ";
  if (configuration->SaveMPQC(filename, mol))
    Log() << Verbose(0) << "done." << endl;
  else
    Log() << Verbose(0) << "failed." << endl;

  if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
  }

  delete(mol);
};

/** Parses the command line options.
 * \param argc argument count
 * \param **argv arguments array
 * \param *molecules list of molecules structure
 * \param *periode elements structure
 * \param configuration config file structure
 * \param *ConfigFileName pointer to config file name in **argv
 * \param *PathToDatabases pointer to db's path in **argv
 * \return exit code (0 - successful, all else - something's wrong)
 */
static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName)
{
  Vector x,y,z,n;  // coordinates for absolute point in cell volume
  double *factor; // unit factor if desired
  ifstream test;
  ofstream output;
  string line;
  atom *first;
  bool SaveFlag = false;
  int ExitFlag = 0;
  int j;
  double volume = 0.;
  enum ConfigStatus configPresent = absent;
  clock_t start,end;
  double MaxDistance = -1;
  int argptr;
  molecule *mol = NULL;
  string BondGraphFileName("\n");
  int verbosity = 0;
  strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);

  if (argc > 1) { // config file specified as option
    // 1. : Parse options that just set variables or print help
    argptr = 1;
    do {
      if (argv[argptr][0] == '-') {
        Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
        argptr++;
        switch(argv[argptr-1][1]) {
          case 'h':
          case 'H':
          case '?':
            Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl;
            Log() << Verbose(0) << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
            Log() << Verbose(0) << "or simply " << argv[0] << " without arguments for interactive session." << endl;
            Log() << Verbose(0) << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
            Log() << Verbose(0) << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
            Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-C <type> [params] <output> <bin output> <BinWidth> <BinStart> <BinEnd>\tPair Correlation analysis." << endl;
            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
            Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
            Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
            Log() << Verbose(0) << "\t-F <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
            Log() << Verbose(0) << "\t-FF <MaxDistance> <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
            Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
            Log() << Verbose(0) << "\t-j\t<path> Store all bonds to file." << endl;
            Log() << Verbose(0) << "\t-J\t<path> Store adjacency per atom to file." << endl;
            Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
            Log() << Verbose(0) << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl;
            Log() << Verbose(0) << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
            Log() << Verbose(0) << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
            Log() << Verbose(0) << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
            Log() << Verbose(0) << "\t-O\tCenter atoms in origin." << endl;
            Log() << Verbose(0) << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
            Log() << Verbose(0) << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
            Log() << Verbose(0) << "\t-r <id>\t\tRemove an atom with given id." << endl;
            Log() << Verbose(0) << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl;
            Log() << Verbose(0) << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-S <file> Store temperatures from the config file in <file>." << endl;
            Log() << Verbose(0) << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
            Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl;
            Log() << Verbose(0) << "\t-V\t\tGives version information." << endl;
            Log() << Verbose(0) << "\t-X\t\tset default name of a molecule." << endl;
            Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl;
            return (1);
            break;
          case 'v':
            while (argv[argptr-1][verbosity+1] == 'v') {
              verbosity++;
            }
            setVerbosity(verbosity);
            Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl;
            break;
          case 'V':
            Log() << Verbose(0) << argv[0] << " " << VERSIONSTRING << endl;
            Log() << Verbose(0) << "Build your own molecule position set." << endl;
            return (1);
            break;
          case 'B':
            if (ExitFlag == 0) ExitFlag = 1;
            if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
              ExitFlag = 255;
              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
              performCriticalExit();
            } else {
              SaveFlag = true;
              j = -1;
              Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
              double * const cell_size = World::get()->cell_size;
              for (int i=0;i<6;i++) {
                cell_size[i] = atof(argv[argptr+i]);
              }
              argptr+=6;
            }
            break;
          case 'e':
            if ((argptr >= argc) || (argv[argptr][0] == '-')) {
              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl);
              performCriticalExit();
            } else {
              Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl;
              strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1);
              argptr+=1;
            }
            break;
          case 'g':
            if ((argptr >= argc) || (argv[argptr][0] == '-')) {
              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl);
              performCriticalExit();
            } else {
              BondGraphFileName = argv[argptr];
              Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl;
              argptr+=1;
            }
            break;
          case 'n':
            Log() << Verbose(0) << "I won't parse trajectories." << endl;
            configuration.FastParsing = true;
            break;
          case 'X':
            {
              char **name = &(World::get()->DefaultName);
              delete[](*name);
              const int length = strlen(argv[argptr]);
              *name = new char[length+2];
              strncpy(*name, argv[argptr], length);
              Log() << Verbose(0) << "Default name of new molecules set to " << *name << "." << endl;
            }
            break;
          default:   // no match? Step on
            argptr++;
            break;
        }
      } else
        argptr++;
    } while (argptr < argc);

    // 3a. Parse the element database
    if (periode->LoadPeriodentafel(configuration.databasepath)) {
      Log() << Verbose(0) << "Element list loaded successfully." << endl;
      //periode->Output();
    } else {
      Log() << Verbose(0) << "Element list loading failed." << endl;
      return 1;
    }
    // 3b. Find config file name and parse if possible, also BondGraphFileName
    if (argv[1][0] != '-') {
      // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
      Log() << Verbose(0) << "Config file given." << endl;
      test.open(argv[1], ios::in);
      if (test == NULL) {
        //return (1);
        output.open(argv[1], ios::out);
        if (output == NULL) {
          Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
          configPresent = absent;
        } else {
          Log() << Verbose(0) << "Empty configuration file." << endl;
          ConfigFileName = argv[1];
          configPresent = empty;
          output.close();
        }
      } else {
        test.close();
        ConfigFileName = argv[1];
        Log() << Verbose(1) << "Specified config file found, parsing ... ";
        switch (configuration.TestSyntax(ConfigFileName, periode)) {
          case 1:
            Log() << Verbose(0) << "new syntax." << endl;
            configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules);
            configPresent = present;
            break;
          case 0:
            Log() << Verbose(0) << "old syntax." << endl;
            configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules);
            configPresent = present;
            break;
          default:
            Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl;
            configPresent = empty;
       }
      }
    } else
      configPresent = absent;
     // set mol to first active molecule
     if (molecules->ListOfMolecules.size() != 0) {
       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
         if ((*ListRunner)->ActiveFlag) {
           mol = *ListRunner;
           break;
         }
     }
     if (mol == NULL) {
       mol = new molecule(periode);
       mol->ActiveFlag = true;
       if (ConfigFileName != NULL)
         mol->SetNameFromFilename(ConfigFileName);
       molecules->insert(mol);
     }
     if (configuration.BG == NULL) {
       configuration.BG = new BondGraph(configuration.GetIsAngstroem());
       if ((!BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) {
         Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
       } else {
         DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
       }
     }

    // 4. parse again through options, now for those depending on elements db and config presence
    argptr = 1;
    do {
      Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl;
      if (argv[argptr][0] == '-') {
        argptr++;
        if ((configPresent == present) || (configPresent == empty)) {
          switch(argv[argptr-1][1]) {
            case 'p':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl;
                if (!mol->AddXYZFile(argv[argptr]))
                  Log() << Verbose(2) << "File not found." << endl;
                else {
                  Log() << Verbose(2) << "File found and parsed." << endl;
                  configPresent = present;
                }
              }
              break;
            case 'a':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3]))) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
                first = new atom;
                first->type = periode->FindElement(atoi(argv[argptr]));
                if (first->type != NULL)
                  Log() << Verbose(2) << "found element " << first->type->name << endl;
                for (int i=NDIM;i--;)
                  first->x.x[i] = atof(argv[argptr+1+i]);
                if (first->type != NULL) {
                  mol->AddAtom(first);  // add to molecule
                  if ((configPresent == empty) && (mol->AtomCount != 0))
                    configPresent = present;
                } else
                  DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the specified element." << endl);
                argptr+=4;
              }
              break;
            default:   // no match? Don't step on (this is done in next switch's default)
              break;
          }
        }
        if (configPresent == present) {
          switch(argv[argptr-1][1]) {
            case 'M':
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl);
                performCriticalExit();
              } else {
                configuration.basis = argv[argptr];
                Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;
                argptr+=1;
              }
              break;
            case 'D':
              if (ExitFlag == 0) ExitFlag = 1;
              {
                Log() << Verbose(1) << "Depth-First-Search Analysis." << endl;
                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
                int *MinimumRingSize = new int[mol->AtomCount];
                atom ***ListOfLocalAtoms = NULL;
                class StackClass<bond *> *BackEdgeStack = NULL;
                class StackClass<bond *> *LocalBackEdgeStack = NULL;
                mol->CreateAdjacencyList(atof(argv[argptr]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
                Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
                if (Subgraphs != NULL) {
                  int FragmentCounter = 0;
                  while (Subgraphs->next != NULL) {
                    Subgraphs = Subgraphs->next;
                    Subgraphs->FillBondStructureFromReference(mol, FragmentCounter, ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
                    LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
                    Subgraphs->Leaf->PickLocalBackEdges(ListOfLocalAtoms[FragmentCounter], BackEdgeStack, LocalBackEdgeStack);
                    Subgraphs->Leaf->CyclicStructureAnalysis(LocalBackEdgeStack, MinimumRingSize);
                    delete(LocalBackEdgeStack);
                    delete(Subgraphs->previous);
                    FragmentCounter++;
                  }
                  delete(Subgraphs);
                  for (int i=0;i<FragmentCounter;i++)
                    Free(&ListOfLocalAtoms[i]);
                  Free(&ListOfLocalAtoms);
                }
                delete(BackEdgeStack);
                delete[](MinimumRingSize);
              }
              //argptr+=1;
              break;
            case 'I':
              Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;
              // @TODO rather do the dissection afterwards
              molecules->DissectMoleculeIntoConnectedSubgraphs(periode, &configuration);
              mol = NULL;
              if (molecules->ListOfMolecules.size() != 0) {
                for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
                  if ((*ListRunner)->ActiveFlag) {
                    mol = *ListRunner;
                    break;
                  }
              }
              if ((mol == NULL) && (!molecules->ListOfMolecules.empty())) {
                mol = *(molecules->ListOfMolecules.begin());
                if (mol != NULL)
                  mol->ActiveFlag = true;
              }
              break;
            case 'C':
              {
                int ranges[3] = {1, 1, 1};
                bool periodic = (argv[argptr-1][2] =='p');
                if (ExitFlag == 0) ExitFlag = 1;
                if ((argptr >= argc)) {
                  ExitFlag = 255;
                  DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C[p] <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl);
                  performCriticalExit();
                } else {
                  switch(argv[argptr][0]) {
                    case 'E':
                      {
                        if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) {
                          ExitFlag = 255;
                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl);
                          performCriticalExit();
                        } else {
                          ofstream output(argv[argptr+3]);
                          ofstream binoutput(argv[argptr+4]);
                          const double BinStart = atof(argv[argptr+5]);
                          const double BinEnd = atof(argv[argptr+6]);

                          element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
                          element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2]));
                          PairCorrelationMap *correlationmap = NULL;
                          if (periodic)
                            correlationmap = PeriodicPairCorrelation(molecules, elemental, elemental2, ranges);
                          else
                            correlationmap = PairCorrelation(molecules, elemental, elemental2);
                          //OutputCorrelationToSurface(&output, correlationmap);
                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
                          OutputCorrelation ( &binoutput, binmap );
                          output.close();
                          binoutput.close();
                          delete(binmap);
                          delete(correlationmap);
                          argptr+=7;
                        }
                      }
                      break;

                    case 'P':
                      {
                        if ((argptr+8 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+7])) || (!IsValidNumber(argv[argptr+8])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) {
                          ExitFlag = 255;
                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl);
                          performCriticalExit();
                        } else {
                          ofstream output(argv[argptr+5]);
                          ofstream binoutput(argv[argptr+6]);
                          const double BinStart = atof(argv[argptr+7]);
                          const double BinEnd = atof(argv[argptr+8]);

                          element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
                          Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3]));
                          CorrelationToPointMap *correlationmap = NULL;
                          if (periodic)
                            correlationmap  = PeriodicCorrelationToPoint(molecules, elemental, Point, ranges);
                          else
                            correlationmap = CorrelationToPoint(molecules, elemental, Point);
                          //OutputCorrelationToSurface(&output, correlationmap);
                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
                          OutputCorrelation ( &binoutput, binmap );
                          output.close();
                          binoutput.close();
                          delete(Point);
                          delete(binmap);
                          delete(correlationmap);
                          argptr+=9;
                        }
                      }
                      break;

                    case 'S':
                      {
                        if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) {
                          ExitFlag = 255;
                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C S <Z> <output> <bin output> <BinWidth> <BinStart> <BinEnd>" << endl);
                          performCriticalExit();
                        } else {
                          ofstream output(argv[argptr+2]);
                          ofstream binoutput(argv[argptr+3]);
                          const double radius = 4.;
                          const double BinWidth = atof(argv[argptr+4]);
                          const double BinStart = atof(argv[argptr+5]);
                          const double BinEnd = atof(argv[argptr+6]);
                          double LCWidth = 20.;
                          if (BinEnd > 0) {
                            if (BinEnd > 2.*radius)
                                LCWidth = BinEnd;
                            else
                              LCWidth = 2.*radius;
                          }

                          // get the boundary
                          class molecule *Boundary = NULL;
                          class Tesselation *TesselStruct = NULL;
                          const LinkedCell *LCList = NULL;
                          // find biggest molecule
                          int counter  = 0;
                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                            if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
                              Boundary = *BigFinder;
                            }
                            counter++;
                          }
                          bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
                          counter = 0;
                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                            Actives[counter++] = (*BigFinder)->ActiveFlag;
                            (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
                          }
                          LCList = new LinkedCell(Boundary, LCWidth);
                          element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
                          FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
                          CorrelationToSurfaceMap *surfacemap = NULL;
                          if (periodic)
                            surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges);
                          else
                            surfacemap = CorrelationToSurface( molecules, elemental, TesselStruct, LCList);
                          OutputCorrelationToSurface(&output, surfacemap);
                          // check whether radius was appropriate
                          {
                            double start; double end;
                            GetMinMax( surfacemap, start, end);
                            if (LCWidth < end)
                              DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl);
                          }
                          BinPairMap *binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd );
                          OutputCorrelation ( &binoutput, binmap );
                          output.close();
                          binoutput.close();
                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
                            (*BigFinder)->ActiveFlag = Actives[counter++];
                          Free(&Actives);
                          delete(LCList);
                          delete(TesselStruct);
                          delete(binmap);
                          delete(surfacemap);
                          argptr+=7;
                        }
                      }
                      break;

                    default:
                      ExitFlag = 255;
                      DoeLog(0) && (eLog()<< Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl);
                      performCriticalExit();
                      break;
                  }
                }
                break;
              }
            case 'E':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
                first = mol->FindAtom(atoi(argv[argptr]));
                first->type = periode->FindElement(atoi(argv[argptr+1]));
                argptr+=2;
              }
              break;
            case 'F':
              if (ExitFlag == 0) ExitFlag = 1;
              MaxDistance = -1;
              if (argv[argptr-1][2] == 'F') { // option is -FF?
                // fetch first argument as max distance to surface
                MaxDistance = atof(argv[argptr++]);
                Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl;
              }
              if ((argptr+7 >=argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+7]))) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <xyz of filler> <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Filling Box with water molecules." << endl;
                // construct water molecule
                molecule *filler = new molecule(periode);
                if (!filler->AddXYZFile(argv[argptr])) {
                  DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse filler molecule from " << argv[argptr] << "." << endl);
                }
                filler->SetNameFromFilename(argv[argptr]);
                configuration.BG->ConstructBondGraph(filler);
                molecule *Filling = NULL;
                // call routine
                double distance[NDIM];
                for (int i=0;i<NDIM;i++)
                  distance[i] = atof(argv[argptr+i+1]);
                Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+4]), atof(argv[argptr+5]), atof(argv[argptr+6]), atoi(argv[argptr+7]));
                if (Filling != NULL) {
                  Filling->ActiveFlag = false;
                  molecules->insert(Filling);
                }
                delete(filler);
                argptr+=6;
              }
              break;
            case 'A':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag =255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl);
                performCriticalExit();
              } else {
                Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl;
                ifstream *input = new ifstream(argv[argptr]);
                mol->CreateAdjacencyListFromDbondFile(input);
                input->close();
                argptr+=1;
              }
              break;

            case 'J':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag =255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of adjacency file: -j <path>" << endl);
                performCriticalExit();
              } else {
                Log() << Verbose(0) << "Storing adjacency to path " << argv[argptr] << "." << endl;
                configuration.BG->ConstructBondGraph(mol);
                mol->StoreAdjacencyToFile(NULL, argv[argptr]);
                argptr+=1;
              }
              break;

            case 'j':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag =255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of bonds file: -j <path>" << endl);
                performCriticalExit();
              } else {
                Log() << Verbose(0) << "Storing bonds to path " << argv[argptr] << "." << endl;
                configuration.BG->ConstructBondGraph(mol);
                mol->StoreBondsToFile(NULL, argv[argptr]);
                argptr+=1;
              }
              break;

            case 'N':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl);
                performCriticalExit();
              } else {
                class Tesselation *T = NULL;
                const LinkedCell *LCList = NULL;
                molecule * Boundary = NULL;
                //string filename(argv[argptr+1]);
                //filename.append(".csv");
                Log() << Verbose(0) << "Evaluating non-convex envelope of biggest molecule.";
                Log() << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
                // find biggest molecule
                int counter  = 0;
                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                  (*BigFinder)->CountAtoms();
                  if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
                    Boundary = *BigFinder;
                  }
                  counter++;
                }
                Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl;
                start = clock();
                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
                if (!FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]))
                  ExitFlag = 255;
                //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
                end = clock();
                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
                delete(LCList);
                delete(T);
                argptr+=2;
              }
              break;
            case 'S':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl);
                performCriticalExit();
              } else {
                Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
                ofstream *output = new ofstream(argv[argptr], ios::trunc);
                if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
                  Log() << Verbose(2) << "File could not be written." << endl;
                else
                  Log() << Verbose(2) << "File stored." << endl;
                output->close();
                delete(output);
                argptr+=1;
              }
              break;
            case 'L':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
                if (atoi(argv[argptr+3]) == 1)
                  Log() << Verbose(1) << "Using Identity for the permutation map." << endl;
                if (!mol->LinearInterpolationBetweenConfiguration(atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false)
                  Log() << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
                else
                  Log() << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
                argptr+=4;
              }
              break;
            case 'P':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
                if (!mol->VerletForceIntegration(argv[argptr], configuration))
                  Log() << Verbose(2) << "File not found." << endl;
                else
                  Log() << Verbose(2) << "File found and parsed." << endl;
                argptr+=1;
              }
              break;
            case 'R':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl;
                double tmp1 = atof(argv[argptr+1]);
                atom *third = mol->FindAtom(atoi(argv[argptr]));
                atom *first = mol->start;
                if ((third != NULL) && (first != mol->end)) {
                  atom *second = first->next;
                  while(second != mol->end) {
                    first = second;
                    second = first->next;
                    if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
                      mol->RemoveAtom(first);
                  }
                } else {
                  DoeLog(1) && (eLog()<< Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl);
                }
                argptr+=2;
              }
              break;
            case 't':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl);
                performCriticalExit();
              } else {
                if (ExitFlag == 0) ExitFlag = 1;
                SaveFlag = true;
                Log() << Verbose(1) << "Translating all ions by given vector." << endl;
                for (int i=NDIM;i--;)
                  x.x[i] = atof(argv[argptr+i]);
                mol->Translate((const Vector *)&x);
                argptr+=3;
              }
              break;
            case 'T':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl);
                performCriticalExit();
              } else {
                if (ExitFlag == 0) ExitFlag = 1;
                SaveFlag = true;
                Log() << Verbose(1) << "Translating all ions periodically by given vector." << endl;
                for (int i=NDIM;i--;)
                  x.x[i] = atof(argv[argptr+i]);
                mol->TranslatePeriodically((const Vector *)&x);
                argptr+=3;
              }
              break;
            case 's':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Scaling all ion positions by factor." << endl;
                factor = new double[NDIM];
                factor[0] = atof(argv[argptr]);
                factor[1] = atof(argv[argptr+1]);
                factor[2] = atof(argv[argptr+2]);
                mol->Scale((const double ** const)&factor);
                double * const cell_size = World::get()->cell_size;
                for (int i=0;i<NDIM;i++) {
                  j += i+1;
                  x.x[i] = atof(argv[NDIM+i]);
                  cell_size[j]*=factor[i];
                }
                delete[](factor);
                argptr+=3;
              }
              break;
            case 'b':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
                double * const cell_size = World::get()->cell_size;
                for (int i=0;i<6;i++) {
                  cell_size[i] = atof(argv[argptr+i]);
                }
                // center
                mol->CenterInBox();
                argptr+=6;
              }
              break;
            case 'B':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
                double * const cell_size = World::get()->cell_size;
                for (int i=0;i<6;i++) {
                  cell_size[i] = atof(argv[argptr+i]);
                }
                // center
                mol->BoundInBox();
                argptr+=6;
              }
              break;
            case 'c':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
                // make every coordinate positive
                mol->CenterEdge(&x);
                // update Box of atoms by boundary
                mol->SetBoxDimension(&x);
                // translate each coordinate by boundary
                double * const cell_size = World::get()->cell_size;
                j=-1;
                for (int i=0;i<NDIM;i++) {
                  j += i+1;
                  x.x[i] = atof(argv[argptr+i]);
                  cell_size[j] += x.x[i]*2.;
                }
                mol->Translate((const Vector *)&x);
                argptr+=3;
              }
              break;
            case 'O':
              if (ExitFlag == 0) ExitFlag = 1;
              SaveFlag = true;
              Log() << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl;
              x.Zero();
              mol->CenterEdge(&x);
              mol->SetBoxDimension(&x);
              argptr+=0;
              break;
            case 'r':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])))  {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl;
                atom *first = mol->FindAtom(atoi(argv[argptr]));
                mol->RemoveAtom(first);
                argptr+=1;
              }
              break;
            case 'f':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl);
                performCriticalExit();
              } else {
                Log() << Verbose(0) << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
                Log() << Verbose(0) << "Creating connection matrix..." << endl;
                start = clock();
                mol->CreateAdjacencyList(atof(argv[argptr++]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
                Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
                if (mol->first->next != mol->last) {
                  ExitFlag = mol->FragmentMolecule(atoi(argv[argptr]), &configuration);
                }
                end = clock();
                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
                argptr+=2;
              }
              break;
            case 'm':
              if (ExitFlag == 0) ExitFlag = 1;
              j = atoi(argv[argptr++]);
              if ((j<0) || (j>1)) {
                DoeLog(1) && (eLog()<< Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl);
                j = 0;
              }
              if (j) {
                SaveFlag = true;
                Log() << Verbose(0) << "Converting to prinicipal axis system." << endl;
              } else
                Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;
              mol->PrincipalAxisSystem((bool)j);
              break;
            case 'o':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl);
                performCriticalExit();
              } else {
                class Tesselation *TesselStruct = NULL;
                const LinkedCell *LCList = NULL;
                Log() << Verbose(0) << "Evaluating volume of the convex envelope.";
                Log() << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;
                Log() << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;
                LCList = new LinkedCell(mol, 10.);
                //FindConvexBorder(mol, LCList, argv[argptr]);
                FindNonConvexBorder(mol, TesselStruct, LCList, 5., argv[argptr+1]);
//                RemoveAllBoundaryPoints(TesselStruct, mol, argv[argptr]);
                double volumedifference = ConvexizeNonconvexEnvelope(TesselStruct, mol, argv[argptr]);
                double clustervolume = VolumeOfConvexEnvelope(TesselStruct, &configuration);
                Log() << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                Log() << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                delete(TesselStruct);
                delete(LCList);
                argptr+=2;
              }
              break;
            case 'U':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl);
                performCriticalExit();
              } else {
                volume = atof(argv[argptr++]);
                Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
              }
            case 'u':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
                if (volume != -1)
                  ExitFlag = 255;
                  DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl);
                  performCriticalExit();
              } else {
                double density;
                SaveFlag = true;
                Log() << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
                density = atof(argv[argptr++]);
                if (density < 1.0) {
                  DoeLog(1) && (eLog()<< Verbose(1) << "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)
//                    DoeLog(1) && (eLog()<< Verbose(1) << "repetition value must be greater 1!" << endl);
//                  repetition[i] = 1;
//                }
                PrepareClustersinWater(&configuration, mol, volume, density);  // if volume == 0, will calculate from ConvexEnvelope
              }
              break;
            case 'd':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl);
                performCriticalExit();
              } else {
                SaveFlag = true;
                double * const cell_size = World::get()->cell_size;
                for (int axis = 1; axis <= NDIM; axis++) {
                  int faktor = atoi(argv[argptr++]);
                  int count;
                  element ** Elements;
                  Vector ** vectors;
                  if (faktor < 1) {
                    DoeLog(1) && (eLog()<< Verbose(1) << "Repetition factor mus be greater than 1!" << endl);
                    faktor = 1;
                  }
                  mol->CountAtoms();  // 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)
                      DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
                    x.Zero();
                    y.Zero();
                    y.x[abs(axis)-1] = 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);
                    }
                    cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
                  }
                }
              }
              break;
            default:   // no match? Step on
              if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
                argptr++;
              break;
          }
        }
      } else argptr++;
    } while (argptr < argc);
    if (SaveFlag)
      SaveConfig(ConfigFileName, &configuration, periode, molecules);
  } else {  // no arguments, hence scan the elements db
    if (periode->LoadPeriodentafel(configuration.databasepath))
      Log() << Verbose(0) << "Element list loaded successfully." << endl;
    else
      Log() << Verbose(0) << "Element list loading failed." << endl;
    configuration.RetrieveConfigPathAndName("main_pcp_linux");
  }
  return(ExitFlag);
};

/********************************************** Main routine **************************************/

int main(int argc, char **argv)
{
  periodentafel *periode = new periodentafel; // and a period table of all elements
  MoleculeListClass *molecules = new MoleculeListClass;  // list of all molecules
  molecule *mol = NULL;
  config *configuration = new config;
  char choice;  // menu choice char
  Vector x,y,z,n;  // coordinates for absolute point in cell volume
  ifstream test;
  ofstream output;
  string line;
  char *ConfigFileName = NULL;
  int j;

  cout << ESPACKVersion << endl;

  Log() << Verbose(1) << "test" << endl;
  DoLog(3) && (Log() << Verbose(1) << "test");

  setVerbosity(0);

  // =========================== PARSE COMMAND LINE OPTIONS ====================================
  j = ParseCommandLineOptions(argc, argv, molecules, periode, *configuration, ConfigFileName);
  switch(j) {
    case 255:  // something went wrong
    case 2:  // just for -f option
    case 1:  // just for -v and -h options
      delete(molecules); // also free's all molecules contained
      delete(periode);
      delete(configuration);
      Log() << Verbose(0) <<  "Maximum of allocated memory: "
        << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
      Log() << Verbose(0) <<  "Remaining non-freed memory: "
        << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
      MemoryUsageObserver::getInstance()->purgeInstance();
      logger::purgeInstance();
      errorLogger::purgeInstance();
     return (j == 1 ? 0 : j);
    default:
      break;
  }

  // General stuff
  if (molecules->ListOfMolecules.size() == 0) {
    mol = new molecule(periode);
    double * const cell_size = World::get()->cell_size;
    if (cell_size[0] == 0.) {
      Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
      for (int i=0;i<6;i++) {
        Log() << Verbose(1) << "Cell size" << i << ": ";
        cin >> cell_size[i];
      }
    }
    mol->ActiveFlag = true;
    molecules->insert(mol);
  }

  // =========================== START INTERACTIVE SESSION ====================================

  // now the main construction loop
  Log() << Verbose(0) << endl << "Now comes the real construction..." << endl;
  do {
    Log() << Verbose(0) << endl << endl;
    Log() << Verbose(0) << "============Molecule list=======================" << endl;
    molecules->Enumerate((ofstream *)&cout);
    Log() << Verbose(0) << "============Menu===============================" << endl;
    Log() << Verbose(0) << "a - set molecule (in)active" << endl;
    Log() << Verbose(0) << "e - edit molecules (load, parse, save)" << endl;
    Log() << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;
    Log() << Verbose(0) << "M - Merge molecules" << endl;
    Log() << Verbose(0) << "m - manipulate atoms" << endl;
    Log() << Verbose(0) << "-----------------------------------------------" << endl;
    Log() << Verbose(0) << "c - edit the current configuration" << endl;
    Log() << Verbose(0) << "-----------------------------------------------" << endl;
    Log() << Verbose(0) << "s - save current setup to config file" << endl;
    Log() << Verbose(0) << "T - call the current test routine" << endl;
    Log() << Verbose(0) << "q - quit" << endl;
    Log() << Verbose(0) << "===============================================" << endl;
    Log() << Verbose(0) << "Input: ";
    cin >> choice;

    switch (choice) {
      case 'a':  // (in)activate molecule
        {
          Log() << Verbose(0) << "Enter index of molecule: ";
          cin >> j;
          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
            if ((*ListRunner)->IndexNr == j)
              (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
        }
        break;

      case 'c': // edit each field of the configuration
       configuration->Edit();
       break;

      case 'e': // create molecule
        EditMolecules(periode, molecules);
        break;

      case 'g': // manipulate molecules
        ManipulateMolecules(periode, molecules, configuration);
        break;

      case 'M':  // merge molecules
        MergeMolecules(periode, molecules);
        break;

      case 'm': // manipulate atoms
        ManipulateAtoms(periode, molecules, configuration);
        break;

      case 'q': // quit
        break;

      case 's': // save to config file
        SaveConfig(ConfigFileName, configuration, periode, molecules);
        break;

      case 'T':
        testroutine(molecules);
        break;

      default:
        break;
    };
  } while (choice != 'q');

  // save element data base
  if (periode->StorePeriodentafel(configuration->databasepath)) //ElementsFileName
    Log() << Verbose(0) << "Saving of elements.db successful." << endl;
  else
    Log() << Verbose(0) << "Saving of elements.db failed." << endl;

  delete(molecules); // also free's all molecules contained
  delete(periode);
  delete(configuration);

  Log() << Verbose(0) <<  "Maximum of allocated memory: "
    << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
  Log() << Verbose(0) <<  "Remaining non-freed memory: "
    << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
  MemoryUsageObserver::purgeInstance();
  logger::purgeInstance();
  errorLogger::purgeInstance();

  return (0);
}

/********************************************** E N D **************************************************/