source: src/Legacy/oldmenu.cpp@ c89fb4

/** \file menu.cpp
 * The class in this file is responsible for displaying the menu and enabling choices.
 *
 * This class is currently being refactored. Functions were copied from builder.cpp and are
 * to be imported into the menu class.
 *
 */

#include "Helpers/MemDebug.hpp"

#include "Legacy/oldmenu.hpp"
#include "analysis_bonds.hpp"
#include "analysis_correlation.hpp"
#include "World.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "bondgraph.hpp"
#include "boundary.hpp"
#include "config.hpp"
#include "element.hpp"
#include "ellipsoid.hpp"
#include "Helpers/helpers.hpp"
#include "leastsquaremin.hpp"
#include "linkedcell.hpp"
#include "Helpers/Log.hpp"
#include "memoryusageobserverunittest.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "vector_ops.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/Line.hpp"

#include "UIElements/UIFactory.hpp"
#include "UIElements/Dialog.hpp"
#include "Menu/Menu.hpp"
#include "Menu/TextMenu.hpp"
#include "Menu/ActionMenuItem.hpp"
#include "Menu/SeperatorItem.hpp"
#include "Menu/DisplayMenuItem.hpp"
#include "Menu/SubMenuItem.hpp"
#include "Actions/MapOfActions.hpp"
#include "Actions/MethodAction.hpp"
#include "Actions/ErrorAction.hpp"
#include "Views/StreamStringView.hpp"
#include "Views/MethodStringView.hpp"


#include <boost/bind.hpp>

/* copied methods for refactoring */
/*TODO: Move these methods inside menu class
 * and restructure menu class*/

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

/** Submenu for adding atoms to the molecule.
 * \param *periode periodentafel
 * \param *molecule molecules with atoms
 */
void oldmenu::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;
  bool aborted;

  Log() << Verbose(0) << "===========ADD ATOM============================" << endl;
  Log() << Verbose(0) << " a - state absolute coordinates of atom" << endl;
  Log() << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
  Log() << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
  Log() << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
  Log() << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
  Log() << Verbose(0) << "all else - go back" << endl;
  Log() << Verbose(0) << "===============================================" << endl;
  Log() << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
  Log() << Verbose(0) << "INPUT: ";
  cin >> choice;

  switch (choice) {
    default:
      eLog() << Verbose(2) << "Not a valid choice." << endl;
      break;
      case 'a': // absolute coordinates of atom
      {
        Dialog *dialog = UIFactory::getInstance().makeDialog();
        first = World::getInstance().createAtom();
        std::vector<element *> elements;
        dialog->queryVector("Please enter coordinates: ",&first->AtomicPosition,World::getInstance().getDomain(), false);
        dialog->queryElement("Please choose element: ",&elements);
        if(dialog->display()){
          if (elements.size() == 1) {
            first->setType(elements.at(0));
            mol->AddAtom(first);  // add to molecule
          } else {
            DoeLog(1) && (eLog() << Verbose(1) << "Unequal to one element given for element of new atom." << endl);
          }
        }
        else{
          World::getInstance().destroyAtom(first);
        }
      }
      break;

      case 'b': // relative coordinates of atom wrt to reference point
        first = World::getInstance().createAtom();
        valid = true;
        aborted = false;
        do {
          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
          auto_ptr<Dialog> dialog(UIFactory::getInstance().makeDialog());
          dialog->queryVector("Enter reference coordinates.",&x,World::getInstance().getDomain(), true);
          dialog->queryVector("Enter relative coordinates.",&first->AtomicPosition,World::getInstance().getDomain(), false);
          if((aborted = !dialog->display())){
            continue;
          }
          first->AtomicPosition += x;
        } while (!aborted && !(valid = mol->CheckBounds((const Vector *)&first->AtomicPosition)));
        if(!aborted){
          first->AtomicElement = periode->AskElement();  // give type
          mol->AddAtom(first);  // add to molecule
        }
        else{
          World::getInstance().destroyAtom(first);
        }
        break;

      case 'c': // relative coordinates of atom wrt to already placed atom
      {
        first = World::getInstance().createAtom();
        valid = true;
        do {
          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
          auto_ptr<Dialog> dialog(UIFactory::getInstance().makeDialog());
          second = mol->AskAtom("Enter atom number: ");
          dialog->queryVector("Enter relative coordinates.",&first->AtomicPosition,World::getInstance().getDomain(), false);
          dialog->display();
          for (int i=NDIM;i--;) {
            first->AtomicPosition[i] += second->AtomicPosition[i];
          }
        } while (!(valid = mol->CheckBounds((const Vector *)&first->AtomicPosition)));
        first->AtomicElement = periode->AskElement();  // give type
        mol->AddAtom(first);  // add to molecule
      }
      break;

    case 'd': // two atoms, two angles and a distance
        first = World::getInstance().createAtom();
        valid = true;
        do {
          if (!valid) {
            eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->AtomicPosition << endl;
          }
          Log() << 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);
          Log() << 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)) {
            Log() << 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= second->AtomicPosition - third->AtomicPosition;
          x.Normalize();
          Log() << Verbose(0) << "x: " << x << endl;
          z = Plane(second->AtomicPosition,third->AtomicPosition,fourth->AtomicPosition).getNormal();
          Log() << Verbose(0) << "z: " << z << endl;
          y = Plane(x,z,0).getNormal();
          Log() << Verbose(0) << "y: " << y << endl;

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

          // rotate another vector around second angle
          n = y;
          n = Line(zeroVec,x).rotateVector(n,c - M_PI);
          Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;

          // add the two linear independent vectors
          first->AtomicPosition += n;
          first->AtomicPosition.Normalize();
          first->AtomicPosition.Scale(a);
          first->AtomicPosition += second->AtomicPosition;

          Log() << Verbose(0) << "resulting coordinates: " << first->AtomicPosition << endl;
        } while (!(valid = mol->CheckBounds((const Vector *)&first->AtomicPosition)));
        first->AtomicElement = periode->AskElement();  // give type
        mol->AddAtom(first);  // add to molecule
        break;

      case 'e': // least square distance position to a set of atoms
        first = World::getInstance().createAtom();
        atoms = new (Vector*[128]);
        valid = true;
        for(int i=128;i--;)
          atoms[i] = NULL;
        int i=0, j=0;
        Log() << Verbose(0) << "Now we need at least three molecules.\n";
        do {
          Log() << Verbose(0) << "Enter " << i+1 << "th atom: ";
          cin >> j;
          if (j != -1) {
            second = mol->FindAtom(j);
            atoms[i++] = &(second->AtomicPosition);
          }
        } while ((j != -1) && (i<128));
        if (i >= 2) {
          LSQdistance(first->AtomicPosition,(const Vector **)atoms, i);

          Log() << Verbose(0) << first->AtomicPosition;
          first->AtomicElement = periode->AskElement();  // give type
          mol->AddAtom(first);  // add to molecule
        } else {
          World::getInstance().destroyAtom(first);
          Log() << 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
 */
void oldmenu::CenterAtoms(molecule *mol)
{
  Vector x, y, helper;
  char choice;  // menu choice char

  Log() << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
  Log() << Verbose(0) << " a - on origin" << endl;
  Log() << Verbose(0) << " b - on center of gravity" << endl;
  Log() << Verbose(0) << " c - within box with additional boundary" << endl;
  Log() << Verbose(0) << " d - within given simulation box" << 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':
      Log() << Verbose(0) << "Centering atoms in config file on origin." << endl;
      mol->CenterOrigin();
      break;
    case 'b':
      Log() << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
      mol->CenterPeriodic();
      break;
    case 'c':
      Log() << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
      for (int i=0;i<NDIM;i++) {
        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
        cin >> y[i];
      }
      mol->CenterEdge(&x);  // make every coordinate positive
      mol->Center += y; // translate by boundary
      helper = (2*y)+x;
      mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
      break;
    case 'd':
      Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
      for (int i=0;i<NDIM;i++) {
        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
        cin >> 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
 */
void oldmenu::AlignAtoms(periodentafel *periode, molecule *mol)
{
  atom *first, *second, *third;
  Vector x,n;
  char choice;  // menu choice char

  Log() << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
  Log() << Verbose(0) << " a - state three atoms defining align plane" << endl;
  Log() << Verbose(0) << " b - state alignment vector" << endl;
  Log() << Verbose(0) << " c - state two atoms in alignment direction" << endl;
  Log() << Verbose(0) << " d - align automatically by least square fit" << 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 = Plane(first->AtomicPosition,second->AtomicPosition,third->AtomicPosition).getNormal();
      break;
    case 'b': // normal vector of mirror plane
    {
      Dialog *dialog = UIFactory::getInstance().makeDialog();
      dialog->queryVector("Enter normal vector of mirror plane.",&n,World::getInstance().getDomain(),false);
      dialog->display();
      delete dialog;
      n.Normalize();
    }
    break;

    case 'c': // three atoms defining mirror plane
      first = mol->AskAtom("Enter first atom: ");
      second = mol->AskAtom("Enter second atom: ");

      n = first->AtomicPosition - second->AtomicPosition;
      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);
      Log() << Verbose(0) << "Element is " << param.type->name << endl;
      mol->GetAlignvector(&param);
      for (int i=NDIM;i--;) {
        x[i] = gsl_vector_get(param.x,i);
        n[i] = gsl_vector_get(param.x,i+NDIM);
      }
      gsl_vector_free(param.x);
      Log() << Verbose(0) << "Offset vector: " << x << endl;
      n.Normalize();
      break;
  };
  Log() << Verbose(0) << "Alignment vector: " << n << endl;
  mol->Align(&n);
};

/** Submenu for mirroring the atoms in the molecule.
 * \param *mol molecule with all the atoms
 */
void oldmenu::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 = Plane(first->AtomicPosition,second->AtomicPosition,third->AtomicPosition).getNormal();
      break;
    case 'b': // normal vector of mirror plane
    {
      Dialog *dialog = UIFactory::getInstance().makeDialog();
      dialog->queryVector("Enter normal vector of mirror plane.",&n,World::getInstance().getDomain(),false);
      dialog->display();
      delete dialog;
      n.Normalize();
    }
    break;

    case 'c': // three atoms defining mirror plane
      first = mol->AskAtom("Enter first atom: ");
      second = mol->AskAtom("Enter second atom: ");

      n = first->AtomicPosition - second->AtomicPosition;
      n.Normalize();
      break;
  };
  Log() << Verbose(0) << "Normal vector: " << n << endl;
  mol->Mirror((const Vector *)&n);
};

/** Submenu for removing the atoms from the molecule.
 * \param *mol molecule with all the atoms
 */
void oldmenu::RemoveAtoms(molecule *mol)
{
  atom *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;
        molecule::iterator runner;
        for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
          runner = iter++;
          if ((*runner)->AtomicPosition.DistanceSquared((*runner)->AtomicPosition) > tmp1*tmp1) // distance to first above radius ...
            mol->RemoveAtom((*runner));
        }
      }
      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;
      molecule::iterator runner;
      for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
        runner = iter++;
        if (((*runner)->AtomicPosition[axis] < tmp1) || ((*runner)->AtomicPosition[axis] > tmp2)) {// out of boundary ...
          //Log() << Verbose(0) << "Atom " << *(*runner) << " with " << (*runner)->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
          mol->RemoveAtom((*runner));
        }
      }
      break;
  };
  //mol->Output();
  choice = 'r';
};

/** Submenu for measuring out the atoms in the molecule.
 * \param *periode periodentafel
 * \param *mol molecule with all the atoms
 */
void oldmenu::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) << " h - count the number of hydrogen bonds" << 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.;

      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
        Z = (*iter)->getType()->Z;
        tmp1 = 0.;
        if (first != (*iter)) {
          x = first->AtomicPosition - (*iter)->AtomicPosition;
          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 " << ((*iter)->nr << ": " << tmp1 << " a.u." << endl;
      }
      for (int i=MAX_ELEMENTS;i--;)
        if (min[i] != 0.) Log() << Verbose(0) << "Minimum Bond length between " << first->getType()->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 = first->AtomicPosition - second->AtomicPosition;
      tmp1 = x.Norm();
      Log() << Verbose(1) << "Distance vector is " << x << "." << "/n"
            << "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 = first->AtomicPosition - second->AtomicPosition;
      y = third->AtomicPosition - second->AtomicPosition;
      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(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.);
        Boundaries *BoundaryPoints = NULL;
        FindConvexBorder(mol, BoundaryPoints 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;
    case 'h':
      {
        int Z1;
        cout << "Please enter first interface element: ";
        cin >> Z1;
        const element * InterfaceElement = World::getInstance().getPeriode()->FindElement(Z1);
        int Z2;
        cout << "Please enter second interface element: ";
        cin >> Z2;
        const element * InterfaceElement2 = World::getInstance().getPeriode()->FindElement(Z2);
        cout << endl << "There are " << CountHydrogenBridgeBonds(World::getInstance().getMolecules(), InterfaceElement, InterfaceElement2) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << " and " << (InterfaceElement2 != 0 ? InterfaceElement2->name : "None") << "." << endl;
      }
      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
 */
void oldmenu::FragmentAtoms(molecule *mol, config *configuration)
{
  int Order1;
  std::string path;
  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;
  DoLog(0) && (Log() << Verbose(0) << "What's the output path and prefix [e.g. /home/foo/BondFragment]: ");
  cin >> path;
  if (mol->hasBondStructure()) {
    start = clock();
    mol->FragmentMolecule(Order1, path);
    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
 */
void oldmenu::ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
{
  atom *first, *second;
  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) << "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)
    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->AtomicPosition[i]-second->AtomicPosition[i])*(first->AtomicPosition[i] - second->AtomicPosition[i]);
        minBond = sqrt(minBond);
        Log() << Verbose(0) << "Current Bond length between " << first->getType()->name << " Atom " << first->nr << " and " << second->getType()->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->AtomicPosition[i] -= (second->AtomicPosition[i]-first->AtomicPosition[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 '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->setType(Z);
        Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->getType()->name << "." << endl;
      }
      break;
  }
};

void oldmenu::duplicateCell(MoleculeListClass *molecules, config *configuration) {
  molecule *mol = NULL;
  int axis,faktor,count,j;
  atom *first = NULL;
  const element **Elements;
  Vector x,y;
  Vector **vectors;

  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;
    if (mol->getAtomCount() != 0) {  // if there is more than none
      count = mol->getAtomCount();  // is changed becausing of adding, thus has to be stored away beforehand
      Elements = new const element *[count];
      vectors = new Vector *[count];
      j = 0;
      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
        Elements[j] = (*iter)->AtomicElement;
        vectors[j] = &(*iter)->AtomicPosition;
        j++;
      }
      if (count != j)
        eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
      x.Zero();
      y.Zero();
      y[abs(axis)-1] = World::getInstance().getDomain()[(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 += y; // per factor one cell width further
        for (int k=count;k--;) { // go through every atom of the original cell
          first = World::getInstance().createAtom(); // create a new body
          first->AtomicPosition = (*vectors[k]) + x;  // use coordinate of original atom
          first->AtomicElement = Elements[k];  // insert original element
          mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
        }
      }
      if (mol->hasBondStructure())
        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 = y;
        x.Scale(-(faktor-1));
        mol->Translate(&x);
      }
      World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    }
  }
}

/** Submenu for manipulating molecules.
 * \param *periode periodentafel
 * \param *molecules list of molecule to manipulate
 */
void oldmenu::ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
{
  Vector x,y,z,n; // coordinates for absolute point in cell volume
  char choice;  // menu choice char
  molecule *mol = NULL;
  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)
    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
      duplicateCell(molecules, configuration);
      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;
        Dialog *dialog = UIFactory::getInstance().makeDialog();
        dialog->queryVector("Enter translation vector.",&x,World::getInstance().getDomain(),false);
        if(dialog->display()){
          mol->Center += x;
        }
        delete dialog;
     }
     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);
  }
};


void oldmenu::SimpleAddMolecules(MoleculeListClass *molecules) {
  molecule *srcmol = NULL, *destmol = NULL;
  Dialog *dialog = UIFactory::getInstance().makeDialog();
  dialog->queryMolecule("molecule-by-id",&destmol, MapOfActions::getInstance().getDescription("molecule-by-id"));
  dialog->queryMolecule("molecule-by-id",&srcmol, MapOfActions::getInstance().getDescription("molecule-by-id"));
  if(dialog->display()) {
    molecules->SimpleAdd(srcmol, destmol);
  }
  else {
    Log() << Verbose(0) << "Adding of molecules canceled" << endl;
  }
  delete dialog;
}

void oldmenu::embeddMolecules(MoleculeListClass *molecules) {
  molecule *srcmol = NULL, *destmol = NULL;
  Dialog *dialog = UIFactory::getInstance().makeDialog();
  dialog->queryMolecule("molecule-by-id",&destmol, MapOfActions::getInstance().getDescription("molecule-by-id"));
  dialog->queryMolecule("molecule-by-id",&srcmol, MapOfActions::getInstance().getDescription("molecule-by-id"));
  if(dialog->display()) {
    molecules->EmbedMerge(destmol, srcmol);
  }
  else {
    Log() << Verbose(0) << "embedding of molecules canceled" << endl;
  }


}

void oldmenu::multiMergeMolecules(MoleculeListClass *molecules) {
  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);
  }
}

void oldmenu::simpleMergeMolecules(MoleculeListClass *molecules) {
  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);
  }
}

/** Submenu for merging molecules.
 * \param *periode periodentafel
 * \param *molecules list of molecules to add to
 */
void oldmenu::MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
{
  TextMenu *MergeMoleculesMenu = new TextMenu(Log() << Verbose(0), "Merge Molecules");

  Action *simpleAddAction = new MethodAction("simpleAddAction",boost::bind(&oldmenu::SimpleAddMolecules,this,molecules),false);
  new ActionMenuItem('a',"simple add of one molecule to another",MergeMoleculesMenu,simpleAddAction);

  Action *embeddAction = new MethodAction("embeddAction",boost::bind(&oldmenu::embeddMolecules,this,molecules),false);
  new ActionMenuItem('e',"embedding merge of two molecules",MergeMoleculesMenu,embeddAction);

  Action *multiMergeAction = new MethodAction("multiMergeAction",boost::bind(&oldmenu::multiMergeMolecules,this,molecules),false);
  new ActionMenuItem('m',"multi-merge of all molecules",MergeMoleculesMenu,multiMergeAction);

  Action *scatterMergeAction = new ErrorAction("scatterMergeAction","Not Implemented yet",false);
  new ActionMenuItem('s',"scatter merge of two molecules",MergeMoleculesMenu,scatterMergeAction);

  Action *simpleMergeAction = new MethodAction("simpleMergeAction",boost::bind(&oldmenu::simpleMergeMolecules,this,molecules),false);
  new ActionMenuItem('t',"simple merge of two molecules",MergeMoleculesMenu,simpleMergeAction);

  Action *returnAction = new MethodAction("returnAction",boost::bind(&TextMenu::doQuit,MergeMoleculesMenu),false);
  MenuItem *returnItem = new ActionMenuItem('q',"return to Main menu",MergeMoleculesMenu,returnAction);

  MergeMoleculesMenu->addDefault(returnItem);

  MergeMoleculesMenu->display();
};


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

/** Is called always as option 'T' in the menu.
 * \param *molecules list of molecules
 */
void oldmenu::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 {
    eLog() << Verbose(0) << "I don't have anything to test on ... ";
    performCriticalExit();
    return;
  }

  // generate some KeySets
  Log() << Verbose(0) << "Generating KeySets." << endl;
  KeySet TestSets[mol->getAtomCount()+1];
  i=1;
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    for (int j=0;j<i;j++) {
      TestSets[j].insert((*iter)->nr);
    }
    i++;
  }
  Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
  KeySetTestPair test;
  molecule::const_iterator iter = mol->begin();
  if (iter != mol->end()) {
    test = TestSets[mol->getAtomCount()-1].insert((*iter)->nr);
    if (test.second) {
      Log() << Verbose(1) << "Insertion worked?!" << endl;
    } else {
      Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
    }
  } else {
    eLog() << Verbose(1) << "No atoms to test double insertion." << endl;
  }

  // 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->getAtomCount();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->getAtomCount()],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++;
  }
  World::getInstance().destroyMolecule(mol);
};

oldmenu::oldmenu()
{
  // TODO Auto-generated constructor stub
}

oldmenu::~oldmenu()
{
  // TODO Auto-generated destructor stub
}