Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    r45f5d6 r244a84  
    4848
    4949
    50 #include <boost/bind.hpp>
    51 
    5250using namespace std;
     51
     52#include <cstring>
    5353
    5454#include "analysis_correlation.hpp"
    … …  
    6767#include "molecule.hpp"
    6868#include "periodentafel.hpp"
    69 #include "UIElements/UIFactory.hpp"
    70 #include "UIElements/MainWindow.hpp"
    71 #include "UIElements/Dialog.hpp"
    72 #include "Menu/ActionMenuItem.hpp"
    73 #include "Actions/ActionRegistry.hpp"
    74 #include "Actions/MethodAction.hpp"
    75 
     69#include "version.h"
     70
     71/********************************************* Subsubmenu routine ************************************/
     72
     73/** Submenu for adding atoms to the molecule.
     74 * \param *periode periodentafel
     75 * \param *molecule molecules with atoms
     76 */
     77static void AddAtoms(periodentafel *periode, molecule *mol)
     78{
     79  atom *first, *second, *third, *fourth;
     80  Vector **atoms;
     81  Vector x,y,z,n;  // coordinates for absolute point in cell volume
     82  double a,b,c;
     83  char choice;  // menu choice char
     84  bool valid;
     85
     86  Log() << Verbose(0) << "===========ADD ATOM============================" << endl;
     87  Log() << Verbose(0) << " a - state absolute coordinates of atom" << endl;
     88  Log() << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
     89  Log() << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
     90  Log() << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
     91  Log() << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
     92  Log() << Verbose(0) << "all else - go back" << endl;
     93  Log() << Verbose(0) << "===============================================" << endl;
     94  Log() << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
     95  Log() << Verbose(0) << "INPUT: ";
     96  cin >> choice;
     97
     98  switch (choice) {
     99    default:
     100      eLog() << Verbose(2) << "Not a valid choice." << endl;
     101      break;
     102      case 'a': // absolute coordinates of atom
     103        Log() << Verbose(0) << "Enter absolute coordinates." << endl;
     104        first = new atom;
     105        first->x.AskPosition(mol->cell_size, false);
     106        first->type = periode->AskElement();  // give type
     107        mol->AddAtom(first);  // add to molecule
     108        break;
     109
     110      case 'b': // relative coordinates of atom wrt to reference point
     111        first = new atom;
     112        valid = true;
     113        do {
     114          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
     115          Log() << Verbose(0) << "Enter reference coordinates." << endl;
     116          x.AskPosition(mol->cell_size, true);
     117          Log() << Verbose(0) << "Enter relative coordinates." << endl;
     118          first->x.AskPosition(mol->cell_size, false);
     119          first->x.AddVector((const Vector *)&x);
     120          Log() << Verbose(0) << "\n";
     121        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     122        first->type = periode->AskElement();  // give type
     123        mol->AddAtom(first);  // add to molecule
     124        break;
     125
     126      case 'c': // relative coordinates of atom wrt to already placed atom
     127        first = new atom;
     128        valid = true;
     129        do {
     130          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
     131          second = mol->AskAtom("Enter atom number: ");
     132          Log() << Verbose(0) << "Enter relative coordinates." << endl;
     133          first->x.AskPosition(mol->cell_size, false);
     134          for (int i=NDIM;i--;) {
     135            first->x.x[i] += second->x.x[i];
     136          }
     137        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     138        first->type = periode->AskElement();  // give type
     139        mol->AddAtom(first);  // add to molecule
     140        break;
     141
     142    case 'd': // two atoms, two angles and a distance
     143        first = new atom;
     144        valid = true;
     145        do {
     146          if (!valid) {
     147            eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl;
     148          }
     149          Log() << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
     150          second = mol->AskAtom("Enter central atom: ");
     151          third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
     152          fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
     153          a = ask_value("Enter distance between central (first) and new atom: ");
     154          b = ask_value("Enter angle between new, first and second atom (degrees): ");
     155          b *= M_PI/180.;
     156          bound(&b, 0., 2.*M_PI);
     157          c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
     158          c *= M_PI/180.;
     159          bound(&c, -M_PI, M_PI);
     160          Log() << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
     161/*
     162          second->Output(1,1,(ofstream *)&cout);
     163          third->Output(1,2,(ofstream *)&cout);
     164          fourth->Output(1,3,(ofstream *)&cout);
     165          n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
     166          x.Copyvector(&second->x);
     167          x.SubtractVector(&third->x);
     168          x.Copyvector(&fourth->x);
     169          x.SubtractVector(&third->x);
     170
     171          if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
     172            Log() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
     173            continue;
     174          }
     175          Log() << Verbose(0) << "resulting relative coordinates: ";
     176          z.Output();
     177          Log() << Verbose(0) << endl;
     178          */
     179          // calc axis vector
     180          x.CopyVector(&second->x);
     181          x.SubtractVector(&third->x);
     182          x.Normalize();
     183          Log() << Verbose(0) << "x: ",
     184          x.Output();
     185          Log() << Verbose(0) << endl;
     186          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
     187          Log() << Verbose(0) << "z: ",
     188          z.Output();
     189          Log() << Verbose(0) << endl;
     190          y.MakeNormalVector(&x,&z);
     191          Log() << Verbose(0) << "y: ",
     192          y.Output();
     193          Log() << Verbose(0) << endl;
     194
     195          // rotate vector around first angle
     196          first->x.CopyVector(&x);
     197          first->x.RotateVector(&z,b - M_PI);
     198          Log() << Verbose(0) << "Rotated vector: ",
     199          first->x.Output();
     200          Log() << Verbose(0) << endl;
     201          // remove the projection onto the rotation plane of the second angle
     202          n.CopyVector(&y);
     203          n.Scale(first->x.ScalarProduct(&y));
     204          Log() << Verbose(0) << "N1: ",
     205          n.Output();
     206          Log() << Verbose(0) << endl;
     207          first->x.SubtractVector(&n);
     208          Log() << Verbose(0) << "Subtracted vector: ",
     209          first->x.Output();
     210          Log() << Verbose(0) << endl;
     211          n.CopyVector(&z);
     212          n.Scale(first->x.ScalarProduct(&z));
     213          Log() << Verbose(0) << "N2: ",
     214          n.Output();
     215          Log() << Verbose(0) << endl;
     216          first->x.SubtractVector(&n);
     217          Log() << Verbose(0) << "2nd subtracted vector: ",
     218          first->x.Output();
     219          Log() << Verbose(0) << endl;
     220
     221          // rotate another vector around second angle
     222          n.CopyVector(&y);
     223          n.RotateVector(&x,c - M_PI);
     224          Log() << Verbose(0) << "2nd Rotated vector: ",
     225          n.Output();
     226          Log() << Verbose(0) << endl;
     227
     228          // add the two linear independent vectors
     229          first->x.AddVector(&n);
     230          first->x.Normalize();
     231          first->x.Scale(a);
     232          first->x.AddVector(&second->x);
     233
     234          Log() << Verbose(0) << "resulting coordinates: ";
     235          first->x.Output();
     236          Log() << Verbose(0) << endl;
     237        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     238        first->type = periode->AskElement();  // give type
     239        mol->AddAtom(first);  // add to molecule
     240        break;
     241
     242      case 'e': // least square distance position to a set of atoms
     243        first = new atom;
     244        atoms = new (Vector*[128]);
     245        valid = true;
     246        for(int i=128;i--;)
     247          atoms[i] = NULL;
     248        int i=0, j=0;
     249        Log() << Verbose(0) << "Now we need at least three molecules.\n";
     250        do {
     251          Log() << Verbose(0) << "Enter " << i+1 << "th atom: ";
     252          cin >> j;
     253          if (j != -1) {
     254            second = mol->FindAtom(j);
     255            atoms[i++] = &(second->x);
     256          }
     257        } while ((j != -1) && (i<128));
     258        if (i >= 2) {
     259          first->x.LSQdistance((const Vector **)atoms, i);
     260
     261          first->x.Output();
     262          first->type = periode->AskElement();  // give type
     263          mol->AddAtom(first);  // add to molecule
     264        } else {
     265          delete first;
     266          Log() << Verbose(0) << "Please enter at least two vectors!\n";
     267        }
     268        break;
     269  };
     270};
     271
     272/** Submenu for centering the atoms in the molecule.
     273 * \param *mol molecule with all the atoms
     274 */
     275static void CenterAtoms(molecule *mol)
     276{
     277  Vector x, y, helper;
     278  char choice;  // menu choice char
     279
     280  Log() << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
     281  Log() << Verbose(0) << " a - on origin" << endl;
     282  Log() << Verbose(0) << " b - on center of gravity" << endl;
     283  Log() << Verbose(0) << " c - within box with additional boundary" << endl;
     284  Log() << Verbose(0) << " d - within given simulation box" << endl;
     285  Log() << Verbose(0) << "all else - go back" << endl;
     286  Log() << Verbose(0) << "===============================================" << endl;
     287  Log() << Verbose(0) << "INPUT: ";
     288  cin >> choice;
     289
     290  switch (choice) {
     291    default:
     292      Log() << Verbose(0) << "Not a valid choice." << endl;
     293      break;
     294    case 'a':
     295      Log() << Verbose(0) << "Centering atoms in config file on origin." << endl;
     296      mol->CenterOrigin();
     297      break;
     298    case 'b':
     299      Log() << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
     300      mol->CenterPeriodic();
     301      break;
     302    case 'c':
     303      Log() << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
     304      for (int i=0;i<NDIM;i++) {
     305        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
     306        cin >> y.x[i];
     307      }
     308      mol->CenterEdge(&x);  // make every coordinate positive
     309      mol->Center.AddVector(&y); // translate by boundary
     310      helper.CopyVector(&y);
     311      helper.Scale(2.);
     312      helper.AddVector(&x);
     313      mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
     314      break;
     315    case 'd':
     316      Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     317      for (int i=0;i<NDIM;i++) {
     318        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
     319        cin >> x.x[i];
     320      }
     321      // update Box of atoms by boundary
     322      mol->SetBoxDimension(&x);
     323      // center
     324      mol->CenterInBox();
     325      break;
     326  }
     327};
     328
     329/** Submenu for aligning the atoms in the molecule.
     330 * \param *periode periodentafel
     331 * \param *mol molecule with all the atoms
     332 */
     333static void AlignAtoms(periodentafel *periode, molecule *mol)
     334{
     335  atom *first, *second, *third;
     336  Vector x,n;
     337  char choice;  // menu choice char
     338
     339  Log() << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
     340  Log() << Verbose(0) << " a - state three atoms defining align plane" << endl;
     341  Log() << Verbose(0) << " b - state alignment vector" << endl;
     342  Log() << Verbose(0) << " c - state two atoms in alignment direction" << endl;
     343  Log() << Verbose(0) << " d - align automatically by least square fit" << endl;
     344  Log() << Verbose(0) << "all else - go back" << endl;
     345  Log() << Verbose(0) << "===============================================" << endl;
     346  Log() << Verbose(0) << "INPUT: ";
     347  cin >> choice;
     348
     349  switch (choice) {
     350    default:
     351    case 'a': // three atoms defining mirror plane
     352      first = mol->AskAtom("Enter first atom: ");
     353      second = mol->AskAtom("Enter second atom: ");
     354      third = mol->AskAtom("Enter third atom: ");
     355
     356      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
     357      break;
     358    case 'b': // normal vector of mirror plane
     359      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
     360      n.AskPosition(mol->cell_size,0);
     361      n.Normalize();
     362      break;
     363    case 'c': // three atoms defining mirror plane
     364      first = mol->AskAtom("Enter first atom: ");
     365      second = mol->AskAtom("Enter second atom: ");
     366
     367      n.CopyVector((const Vector *)&first->x);
     368      n.SubtractVector((const Vector *)&second->x);
     369      n.Normalize();
     370      break;
     371    case 'd':
     372      char shorthand[4];
     373      Vector a;
     374      struct lsq_params param;
     375      do {
     376        fprintf(stdout, "Enter the element of atoms to be chosen: ");
     377        fscanf(stdin, "%3s", shorthand);
     378      } while ((param.type = periode->FindElement(shorthand)) == NULL);
     379      Log() << Verbose(0) << "Element is " << param.type->name << endl;
     380      mol->GetAlignvector(&param);
     381      for (int i=NDIM;i--;) {
     382        x.x[i] = gsl_vector_get(param.x,i);
     383        n.x[i] = gsl_vector_get(param.x,i+NDIM);
     384      }
     385      gsl_vector_free(param.x);
     386      Log() << Verbose(0) << "Offset vector: ";
     387      x.Output();
     388      Log() << Verbose(0) << endl;
     389      n.Normalize();
     390      break;
     391  };
     392  Log() << Verbose(0) << "Alignment vector: ";
     393  n.Output();
     394  Log() << Verbose(0) << endl;
     395  mol->Align(&n);
     396};
     397
     398/** Submenu for mirroring the atoms in the molecule.
     399 * \param *mol molecule with all the atoms
     400 */
     401static void MirrorAtoms(molecule *mol)
     402{
     403  atom *first, *second, *third;
     404  Vector n;
     405  char choice;  // menu choice char
     406
     407  Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
     408  Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
     409  Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl;
     410  Log() << Verbose(0) << " c - state two atoms in normal direction" << endl;
     411  Log() << Verbose(0) << "all else - go back" << endl;
     412  Log() << Verbose(0) << "===============================================" << endl;
     413  Log() << Verbose(0) << "INPUT: ";
     414  cin >> choice;
     415
     416  switch (choice) {
     417    default:
     418    case 'a': // three atoms defining mirror plane
     419      first = mol->AskAtom("Enter first atom: ");
     420      second = mol->AskAtom("Enter second atom: ");
     421      third = mol->AskAtom("Enter third atom: ");
     422
     423      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
     424      break;
     425    case 'b': // normal vector of mirror plane
     426      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
     427      n.AskPosition(mol->cell_size,0);
     428      n.Normalize();
     429      break;
     430    case 'c': // three atoms defining mirror plane
     431      first = mol->AskAtom("Enter first atom: ");
     432      second = mol->AskAtom("Enter second atom: ");
     433
     434      n.CopyVector((const Vector *)&first->x);
     435      n.SubtractVector((const Vector *)&second->x);
     436      n.Normalize();
     437      break;
     438  };
     439  Log() << Verbose(0) << "Normal vector: ";
     440  n.Output();
     441  Log() << Verbose(0) << endl;
     442  mol->Mirror((const Vector *)&n);
     443};
     444
     445/** Submenu for removing the atoms from the molecule.
     446 * \param *mol molecule with all the atoms
     447 */
     448static void RemoveAtoms(molecule *mol)
     449{
     450  atom *first, *second;
     451  int axis;
     452  double tmp1, tmp2;
     453  char choice;  // menu choice char
     454
     455  Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
     456  Log() << Verbose(0) << " a - state atom for removal by number" << endl;
     457  Log() << Verbose(0) << " b - keep only in radius around atom" << endl;
     458  Log() << Verbose(0) << " c - remove this with one axis greater value" << endl;
     459  Log() << Verbose(0) << "all else - go back" << endl;
     460  Log() << Verbose(0) << "===============================================" << endl;
     461  Log() << Verbose(0) << "INPUT: ";
     462  cin >> choice;
     463
     464  switch (choice) {
     465    default:
     466    case 'a':
     467      if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
     468        Log() << Verbose(1) << "Atom removed." << endl;
     469      else
     470        Log() << Verbose(1) << "Atom not found." << endl;
     471      break;
     472    case 'b':
     473      second = mol->AskAtom("Enter number of atom as reference point: ");
     474      Log() << Verbose(0) << "Enter radius: ";
     475      cin >> tmp1;
     476      first = mol->start;
     477      second = first->next;
     478      while(second != mol->end) {
     479        first = second;
     480        second = first->next;
     481        if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
     482          mol->RemoveAtom(first);
     483      }
     484      break;
     485    case 'c':
     486      Log() << Verbose(0) << "Which axis is it: ";
     487      cin >> axis;
     488      Log() << Verbose(0) << "Lower boundary: ";
     489      cin >> tmp1;
     490      Log() << Verbose(0) << "Upper boundary: ";
     491      cin >> tmp2;
     492      first = mol->start;
     493      second = first->next;
     494      while(second != mol->end) {
     495        first = second;
     496        second = first->next;
     497        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
     498          //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
     499          mol->RemoveAtom(first);
     500        }
     501      }
     502      break;
     503  };
     504  //mol->Output();
     505  choice = 'r';
     506};
     507
     508/** Submenu for measuring out the atoms in the molecule.
     509 * \param *periode periodentafel
     510 * \param *mol molecule with all the atoms
     511 */
     512static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
     513{
     514  atom *first, *second, *third;
     515  Vector x,y;
     516  double min[256], tmp1, tmp2, tmp3;
     517  int Z;
     518  char choice;  // menu choice char
     519
     520  Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
     521  Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
     522  Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl;
     523  Log() << Verbose(0) << " c - calculate bond angle" << endl;
     524  Log() << Verbose(0) << " d - calculate principal axis of the system" << endl;
     525  Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
     526  Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;
     527  Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
     528  Log() << Verbose(0) << "all else - go back" << endl;
     529  Log() << Verbose(0) << "===============================================" << endl;
     530  Log() << Verbose(0) << "INPUT: ";
     531  cin >> choice;
     532
     533  switch(choice) {
     534    default:
     535      Log() << Verbose(1) << "Not a valid choice." << endl;
     536      break;
     537    case 'a':
     538      first = mol->AskAtom("Enter first atom: ");
     539      for (int i=MAX_ELEMENTS;i--;)
     540        min[i] = 0.;
     541
     542      second = mol->start;
     543      while ((second->next != mol->end)) {
     544        second = second->next; // advance
     545        Z = second->type->Z;
     546        tmp1 = 0.;
     547        if (first != second) {
     548          x.CopyVector((const Vector *)&first->x);
     549          x.SubtractVector((const Vector *)&second->x);
     550          tmp1 = x.Norm();
     551        }
     552        if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
     553        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
     554      }
     555      for (int i=MAX_ELEMENTS;i--;)
     556        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;
     557      break;
     558
     559    case 'b':
     560      first = mol->AskAtom("Enter first atom: ");
     561      second = mol->AskAtom("Enter second atom: ");
     562      for (int i=NDIM;i--;)
     563        min[i] = 0.;
     564      x.CopyVector((const Vector *)&first->x);
     565      x.SubtractVector((const Vector *)&second->x);
     566      tmp1 = x.Norm();
     567      Log() << Verbose(1) << "Distance vector is ";
     568      x.Output();
     569      Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
     570      break;
     571
     572    case 'c':
     573      Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
     574      first = mol->AskAtom("Enter first atom: ");
     575      second = mol->AskAtom("Enter central atom: ");
     576      third  = mol->AskAtom("Enter last atom: ");
     577      tmp1 = tmp2 = tmp3 = 0.;
     578      x.CopyVector((const Vector *)&first->x);
     579      x.SubtractVector((const Vector *)&second->x);
     580      y.CopyVector((const Vector *)&third->x);
     581      y.SubtractVector((const Vector *)&second->x);
     582      Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     583      Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
     584      break;
     585    case 'd':
     586      Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;
     587      Log() << Verbose(0) << "Shall we rotate? [0/1]: ";
     588      cin >> Z;
     589      if ((Z >=0) && (Z <=1))
     590        mol->PrincipalAxisSystem((bool)Z);
     591      else
     592        mol->PrincipalAxisSystem(false);
     593      break;
     594    case 'e':
     595      {
     596        Log() << Verbose(0) << "Evaluating volume of the convex envelope.";
     597        class Tesselation *TesselStruct = NULL;
     598        const LinkedCell *LCList = NULL;
     599        LCList = new LinkedCell(mol, 10.);
     600        FindConvexBorder(mol, TesselStruct, LCList, NULL);
     601        double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
     602        Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;\
     603        delete(LCList);
     604        delete(TesselStruct);
     605      }
     606      break;
     607    case 'f':
     608      mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps);
     609      break;
     610    case 'g':
     611      {
     612        char filename[255];
     613        Log() << Verbose(0) << "Please enter filename: " << endl;
     614        cin >> filename;
     615        Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
     616        ofstream *output = new ofstream(filename, ios::trunc);
     617        if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
     618          Log() << Verbose(2) << "File could not be written." << endl;
     619        else
     620          Log() << Verbose(2) << "File stored." << endl;
     621        output->close();
     622        delete(output);
     623      }
     624      break;
     625  }
     626};
     627
     628/** Submenu for measuring out the atoms in the molecule.
     629 * \param *mol molecule with all the atoms
     630 * \param *configuration configuration structure for the to be written config files of all fragments
     631 */
     632static void FragmentAtoms(molecule *mol, config *configuration)
     633{
     634  int Order1;
     635  clock_t start, end;
     636
     637  Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
     638  Log() << Verbose(0) << "What's the desired bond order: ";
     639  cin >> Order1;
     640  if (mol->first->next != mol->last) {  // there are bonds
     641    start = clock();
     642    mol->FragmentMolecule(Order1, configuration);
     643    end = clock();
     644    Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     645  } else
     646    Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
     647};
     648
     649/********************************************** Submenu routine **************************************/
     650
     651/** Submenu for manipulating atoms.
     652 * \param *periode periodentafel
     653 * \param *molecules list of molecules whose atoms are to be manipulated
     654 */
     655static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
     656{
     657  atom *first, *second;
     658  molecule *mol = NULL;
     659  Vector x,y,z,n; // coordinates for absolute point in cell volume
     660  double *factor; // unit factor if desired
     661  double bond, minBond;
     662  char choice;  // menu choice char
     663  bool valid;
     664
     665  Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl;
     666  Log() << Verbose(0) << "a - add an atom" << endl;
     667  Log() << Verbose(0) << "r - remove an atom" << endl;
     668  Log() << Verbose(0) << "b - scale a bond between atoms" << endl;
     669  Log() << Verbose(0) << "u - change an atoms element" << endl;
     670  Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
     671  Log() << Verbose(0) << "all else - go back" << endl;
     672  Log() << Verbose(0) << "===============================================" << endl;
     673  if (molecules->NumberOfActiveMolecules() > 1)
     674    eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;
     675  Log() << Verbose(0) << "INPUT: ";
     676  cin >> choice;
     677
     678  switch (choice) {
     679    default:
     680      Log() << Verbose(0) << "Not a valid choice." << endl;
     681      break;
     682
     683    case 'a': // add atom
     684      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     685        if ((*ListRunner)->ActiveFlag) {
     686        mol = *ListRunner;
     687        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     688        AddAtoms(periode, mol);
     689      }
     690      break;
     691
     692    case 'b': // scale a bond
     693      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     694        if ((*ListRunner)->ActiveFlag) {
     695        mol = *ListRunner;
     696        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     697        Log() << Verbose(0) << "Scaling bond length between two atoms." << endl;
     698        first = mol->AskAtom("Enter first (fixed) atom: ");
     699        second = mol->AskAtom("Enter second (shifting) atom: ");
     700        minBond = 0.;
     701        for (int i=NDIM;i--;)
     702          minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
     703        minBond = sqrt(minBond);
     704        Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;
     705        Log() << Verbose(0) << "Enter new bond length [a.u.]: ";
     706        cin >> bond;
     707        for (int i=NDIM;i--;) {
     708          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
     709        }
     710        //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
     711        //second->Output(second->type->No, 1);
     712      }
     713      break;
     714
     715    case 'c': // unit scaling of the metric
     716      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     717        if ((*ListRunner)->ActiveFlag) {
     718        mol = *ListRunner;
     719        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     720       Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
     721       Log() << Verbose(0) << "Enter three factors: ";
     722       factor = new double[NDIM];
     723       cin >> factor[0];
     724       cin >> factor[1];
     725       cin >> factor[2];
     726       valid = true;
     727       mol->Scale((const double ** const)&factor);
     728       delete[](factor);
     729      }
     730     break;
     731
     732    case 'l': // measure distances or angles
     733      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     734        if ((*ListRunner)->ActiveFlag) {
     735        mol = *ListRunner;
     736        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     737        MeasureAtoms(periode, mol, configuration);
     738      }
     739      break;
     740
     741    case 'r': // remove atom
     742      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     743        if ((*ListRunner)->ActiveFlag) {
     744        mol = *ListRunner;
     745        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     746        RemoveAtoms(mol);
     747      }
     748      break;
     749
     750    case 'u': // change an atom's element
     751      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     752        if ((*ListRunner)->ActiveFlag) {
     753        int Z;
     754        mol = *ListRunner;
     755        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     756        first = NULL;
     757        do {
     758          Log() << Verbose(0) << "Change the element of which atom: ";
     759          cin >> Z;
     760        } while ((first = mol->FindAtom(Z)) == NULL);
     761        Log() << Verbose(0) << "New element by atomic number Z: ";
     762        cin >> Z;
     763        first->type = periode->FindElement(Z);
     764        Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
     765      }
     766      break;
     767  }
     768};
     769
     770/** Submenu for manipulating molecules.
     771 * \param *periode periodentafel
     772 * \param *molecules list of molecule to manipulate
     773 */
     774static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
     775{
     776  atom *first = NULL;
     777  Vector x,y,z,n; // coordinates for absolute point in cell volume
     778  int j, axis, count, faktor;
     779  char choice;  // menu choice char
     780  molecule *mol = NULL;
     781  element **Elements;
     782  Vector **vectors;
     783  MoleculeLeafClass *Subgraphs = NULL;
     784
     785  Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;
     786  Log() << Verbose(0) << "c - scale by unit transformation" << endl;
     787  Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
     788  Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
     789  Log() << Verbose(0) << "g - center atoms in box" << endl;
     790  Log() << Verbose(0) << "i - realign molecule" << endl;
     791  Log() << Verbose(0) << "m - mirror all molecules" << endl;
     792  Log() << Verbose(0) << "o - create connection matrix" << endl;
     793  Log() << Verbose(0) << "t - translate molecule by vector" << endl;
     794  Log() << Verbose(0) << "all else - go back" << endl;
     795  Log() << Verbose(0) << "===============================================" << endl;
     796  if (molecules->NumberOfActiveMolecules() > 1)
     797    eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;
     798  Log() << Verbose(0) << "INPUT: ";
     799  cin >> choice;
     800
     801  switch (choice) {
     802    default:
     803      Log() << Verbose(0) << "Not a valid choice." << endl;
     804      break;
     805
     806    case 'd': // duplicate the periodic cell along a given axis, given times
     807      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     808        if ((*ListRunner)->ActiveFlag) {
     809        mol = *ListRunner;
     810        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     811        Log() << Verbose(0) << "State the axis [(+-)123]: ";
     812        cin >> axis;
     813        Log() << Verbose(0) << "State the factor: ";
     814        cin >> faktor;
     815
     816        mol->CountAtoms(); // recount atoms
     817        if (mol->AtomCount != 0) {  // if there is more than none
     818          count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
     819          Elements = new element *[count];
     820          vectors = new Vector *[count];
     821          j = 0;
     822          first = mol->start;
     823          while (first->next != mol->end) { // make a list of all atoms with coordinates and element
     824            first = first->next;
     825            Elements[j] = first->type;
     826            vectors[j] = &first->x;
     827            j++;
     828          }
     829          if (count != j)
     830            eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
     831          x.Zero();
     832          y.Zero();
     833          y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     834          for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
     835            x.AddVector(&y); // per factor one cell width further
     836            for (int k=count;k--;) { // go through every atom of the original cell
     837              first = new atom(); // create a new body
     838              first->x.CopyVector(vectors[k]);  // use coordinate of original atom
     839              first->x.AddVector(&x);     // translate the coordinates
     840              first->type = Elements[k];  // insert original element
     841              mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
     842            }
     843          }
     844          if (mol->first->next != mol->last) // if connect matrix is present already, redo it
     845            mol->CreateAdjacencyList(mol->BondDistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
     846          // free memory
     847          delete[](Elements);
     848          delete[](vectors);
     849          // correct cell size
     850          if (axis < 0) { // if sign was negative, we have to translate everything
     851            x.Zero();
     852            x.AddVector(&y);
     853            x.Scale(-(faktor-1));
     854            mol->Translate(&x);
     855          }
     856          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     857        }
     858      }
     859      break;
     860
     861    case 'f':
     862      FragmentAtoms(mol, configuration);
     863      break;
     864
     865    case 'g': // center the atoms
     866      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     867        if ((*ListRunner)->ActiveFlag) {
     868        mol = *ListRunner;
     869        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     870        CenterAtoms(mol);
     871      }
     872      break;
     873
     874    case 'i': // align all atoms
     875      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     876        if ((*ListRunner)->ActiveFlag) {
     877        mol = *ListRunner;
     878        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     879        AlignAtoms(periode, mol);
     880      }
     881      break;
     882
     883    case 'm': // mirror atoms along a given axis
     884      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     885        if ((*ListRunner)->ActiveFlag) {
     886        mol = *ListRunner;
     887        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     888        MirrorAtoms(mol);
     889      }
     890      break;
     891
     892    case 'o': // create the connection matrix
     893      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     894        if ((*ListRunner)->ActiveFlag) {
     895          mol = *ListRunner;
     896          double bonddistance;
     897          clock_t start,end;
     898          Log() << Verbose(0) << "What's the maximum bond distance: ";
     899          cin >> bonddistance;
     900          start = clock();
     901          mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
     902          end = clock();
     903          Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     904        }
     905      break;
     906
     907    case 't': // translate all atoms
     908      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     909        if ((*ListRunner)->ActiveFlag) {
     910        mol = *ListRunner;
     911        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     912        Log() << Verbose(0) << "Enter translation vector." << endl;
     913        x.AskPosition(mol->cell_size,0);
     914        mol->Center.AddVector((const Vector *)&x);
     915     }
     916     break;
     917  }
     918  // Free all
     919  if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
     920    while (Subgraphs->next != NULL) {
     921      Subgraphs = Subgraphs->next;
     922      delete(Subgraphs->previous);
     923    }
     924    delete(Subgraphs);
     925  }
     926};
     927
     928
     929/** Submenu for creating new molecules.
     930 * \param *periode periodentafel
     931 * \param *molecules list of molecules to add to
     932 */
     933static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
     934{
     935  char choice;  // menu choice char
     936  Vector center;
     937  int nr, count;
     938  molecule *mol = NULL;
     939
     940  Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl;
     941  Log() << Verbose(0) << "c - create new molecule" << endl;
     942  Log() << Verbose(0) << "l - load molecule from xyz file" << endl;
     943  Log() << Verbose(0) << "n - change molecule's name" << endl;
     944  Log() << Verbose(0) << "N - give molecules filename" << endl;
     945  Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl;
     946  Log() << Verbose(0) << "r - remove a molecule" << endl;
     947  Log() << Verbose(0) << "all else - go back" << endl;
     948  Log() << Verbose(0) << "===============================================" << endl;
     949  Log() << Verbose(0) << "INPUT: ";
     950  cin >> choice;
     951
     952  switch (choice) {
     953    default:
     954      Log() << Verbose(0) << "Not a valid choice." << endl;
     955      break;
     956    case 'c':
     957      mol = new molecule(periode);
     958      molecules->insert(mol);
     959      break;
     960
     961    case 'l': // load from XYZ file
     962      {
     963        char filename[MAXSTRINGSIZE];
     964        Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     965        mol = new molecule(periode);
     966        do {
     967          Log() << Verbose(0) << "Enter file name: ";
     968          cin >> filename;
     969        } while (!mol->AddXYZFile(filename));
     970        mol->SetNameFromFilename(filename);
     971        // center at set box dimensions
     972        mol->CenterEdge(&center);
     973        mol->cell_size[0] = center.x[0];
     974        mol->cell_size[1] = 0;
     975        mol->cell_size[2] = center.x[1];
     976        mol->cell_size[3] = 0;
     977        mol->cell_size[4] = 0;
     978        mol->cell_size[5] = center.x[2];
     979        molecules->insert(mol);
     980      }
     981      break;
     982
     983    case 'n':
     984      {
     985        char filename[MAXSTRINGSIZE];
     986        do {
     987          Log() << Verbose(0) << "Enter index of molecule: ";
     988          cin >> nr;
     989          mol = molecules->ReturnIndex(nr);
     990        } while (mol == NULL);
     991        Log() << Verbose(0) << "Enter name: ";
     992        cin >> filename;
     993        strcpy(mol->name, filename);
     994      }
     995      break;
     996
     997    case 'N':
     998      {
     999        char filename[MAXSTRINGSIZE];
     1000        do {
     1001          Log() << Verbose(0) << "Enter index of molecule: ";
     1002          cin >> nr;
     1003          mol = molecules->ReturnIndex(nr);
     1004        } while (mol == NULL);
     1005        Log() << Verbose(0) << "Enter name: ";
     1006        cin >> filename;
     1007        mol->SetNameFromFilename(filename);
     1008      }
     1009      break;
     1010
     1011    case 'p': // parse XYZ file
     1012      {
     1013        char filename[MAXSTRINGSIZE];
     1014        mol = NULL;
     1015        do {
     1016          Log() << Verbose(0) << "Enter index of molecule: ";
     1017          cin >> nr;
     1018          mol = molecules->ReturnIndex(nr);
     1019        } while (mol == NULL);
     1020        Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     1021        do {
     1022          Log() << Verbose(0) << "Enter file name: ";
     1023          cin >> filename;
     1024        } while (!mol->AddXYZFile(filename));
     1025        mol->SetNameFromFilename(filename);
     1026      }
     1027      break;
     1028
     1029    case 'r':
     1030      Log() << Verbose(0) << "Enter index of molecule: ";
     1031      cin >> nr;
     1032      count = 1;
     1033      for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1034        if (nr == (*ListRunner)->IndexNr) {
     1035          mol = *ListRunner;
     1036          molecules->ListOfMolecules.erase(ListRunner);
     1037          delete(mol);
     1038          break;
     1039        }
     1040      break;
     1041  }
     1042};
     1043
     1044
     1045/** Submenu for merging molecules.
     1046 * \param *periode periodentafel
     1047 * \param *molecules list of molecules to add to
     1048 */
     1049static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
     1050{
     1051  char choice;  // menu choice char
     1052
     1053  Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl;
     1054  Log() << Verbose(0) << "a - simple add of one molecule to another" << endl;
     1055  Log() << Verbose(0) << "e - embedding merge of two molecules" << endl;
     1056  Log() << Verbose(0) << "m - multi-merge of all molecules" << endl;
     1057  Log() << Verbose(0) << "s - scatter merge of two molecules" << endl;
     1058  Log() << Verbose(0) << "t - simple merge of two molecules" << endl;
     1059  Log() << Verbose(0) << "all else - go back" << endl;
     1060  Log() << Verbose(0) << "===============================================" << endl;
     1061  Log() << Verbose(0) << "INPUT: ";
     1062  cin >> choice;
     1063
     1064  switch (choice) {
     1065    default:
     1066      Log() << Verbose(0) << "Not a valid choice." << endl;
     1067      break;
     1068
     1069    case 'a':
     1070      {
     1071        int src, dest;
     1072        molecule *srcmol = NULL, *destmol = NULL;
     1073        {
     1074          do {
     1075            Log() << Verbose(0) << "Enter index of destination molecule: ";
     1076            cin >> dest;
     1077            destmol = molecules->ReturnIndex(dest);
     1078          } while ((destmol == NULL) && (dest != -1));
     1079          do {
     1080            Log() << Verbose(0) << "Enter index of source molecule to add from: ";
     1081            cin >> src;
     1082            srcmol = molecules->ReturnIndex(src);
     1083          } while ((srcmol == NULL) && (src != -1));
     1084          if ((src != -1) && (dest != -1))
     1085            molecules->SimpleAdd(srcmol, destmol);
     1086        }
     1087      }
     1088      break;
     1089
     1090    case 'e':
     1091      {
     1092        int src, dest;
     1093        molecule *srcmol = NULL, *destmol = NULL;
     1094        do {
     1095          Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ";
     1096          cin >> src;
     1097          srcmol = molecules->ReturnIndex(src);
     1098        } while ((srcmol == NULL) && (src != -1));
     1099        do {
     1100          Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ";
     1101          cin >> dest;
     1102          destmol = molecules->ReturnIndex(dest);
     1103        } while ((destmol == NULL) && (dest != -1));
     1104        if ((src != -1) && (dest != -1))
     1105          molecules->EmbedMerge(destmol, srcmol);
     1106      }
     1107      break;
     1108
     1109    case 'm':
     1110      {
     1111        int nr;
     1112        molecule *mol = NULL;
     1113        do {
     1114          Log() << Verbose(0) << "Enter index of molecule to merge into: ";
     1115          cin >> nr;
     1116          mol = molecules->ReturnIndex(nr);
     1117        } while ((mol == NULL) && (nr != -1));
     1118        if (nr != -1) {
     1119          int N = molecules->ListOfMolecules.size()-1;
     1120          int *src = new int(N);
     1121          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1122            if ((*ListRunner)->IndexNr != nr)
     1123              src[N++] = (*ListRunner)->IndexNr;       
     1124          molecules->SimpleMultiMerge(mol, src, N);
     1125          delete[](src);
     1126        }
     1127      }
     1128      break;
     1129
     1130    case 's':
     1131      Log() << Verbose(0) << "Not implemented yet." << endl;
     1132      break;
     1133
     1134    case 't':
     1135      {
     1136        int src, dest;
     1137        molecule *srcmol = NULL, *destmol = NULL;
     1138        {
     1139          do {
     1140            Log() << Verbose(0) << "Enter index of destination molecule: ";
     1141            cin >> dest;
     1142            destmol = molecules->ReturnIndex(dest);
     1143          } while ((destmol == NULL) && (dest != -1));
     1144          do {
     1145            Log() << Verbose(0) << "Enter index of source molecule to merge into: ";
     1146            cin >> src;
     1147            srcmol = molecules->ReturnIndex(src);
     1148          } while ((srcmol == NULL) && (src != -1));
     1149          if ((src != -1) && (dest != -1))
     1150            molecules->SimpleMerge(srcmol, destmol);
     1151        }
     1152      }
     1153      break;
     1154  }
     1155};
     1156
     1157
     1158/********************************************** Test routine **************************************/
     1159
     1160/** Is called always as option 'T' in the menu.
     1161 * \param *molecules list of molecules
     1162 */
     1163static void testroutine(MoleculeListClass *molecules)
     1164{
     1165  // the current test routine checks the functionality of the KeySet&Graph concept:
     1166  // We want to have a multiindex (the KeySet) describing a unique subgraph
     1167  int i, comp, counter=0;
     1168
     1169  // create a clone
     1170  molecule *mol = NULL;
     1171  if (molecules->ListOfMolecules.size() != 0) // clone
     1172    mol = (molecules->ListOfMolecules.front())->CopyMolecule();
     1173  else {
     1174    eLog() << Verbose(0) << "I don't have anything to test on ... ";
     1175    performCriticalExit();
     1176    return;
     1177  }
     1178  atom *Walker = mol->start;
     1179
     1180  // generate some KeySets
     1181  Log() << Verbose(0) << "Generating KeySets." << endl;
     1182  KeySet TestSets[mol->AtomCount+1];
     1183  i=1;
     1184  while (Walker->next != mol->end) {
     1185    Walker = Walker->next;
     1186    for (int j=0;j<i;j++) {
     1187      TestSets[j].insert(Walker->nr);
     1188    }
     1189    i++;
     1190  }
     1191  Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
     1192  KeySetTestPair test;
     1193  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
     1194  if (test.second) {
     1195    Log() << Verbose(1) << "Insertion worked?!" << endl;
     1196  } else {
     1197    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
     1198  }
     1199  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
     1200  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
     1201
     1202  // constructing Graph structure
     1203  Log() << Verbose(0) << "Generating Subgraph class." << endl;
     1204  Graph Subgraphs;
     1205
     1206  // insert KeySets into Subgraphs
     1207  Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl;
     1208  for (int j=0;j<mol->AtomCount;j++) {
     1209    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
     1210  }
     1211  Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl;
     1212  GraphTestPair test2;
     1213  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
     1214  if (test2.second) {
     1215    Log() << Verbose(1) << "Insertion worked?!" << endl;
     1216  } else {
     1217    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
     1218  }
     1219
     1220  // show graphs
     1221  Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl;
     1222  Graph::iterator A = Subgraphs.begin();
     1223  while (A !=  Subgraphs.end()) {
     1224    Log() << Verbose(0) << (*A).second.first << ": ";
     1225    KeySet::iterator key = (*A).first.begin();
     1226    comp = -1;
     1227    while (key != (*A).first.end()) {
     1228      if ((*key) > comp)
     1229        Log() << Verbose(0) << (*key) << " ";
     1230      else
     1231        Log() << Verbose(0) << (*key) << "! ";
     1232      comp = (*key);
     1233      key++;
     1234    }
     1235    Log() << Verbose(0) << endl;
     1236    A++;
     1237  }
     1238  delete(mol);
     1239};
     1240
     1241/** Tries given filename or standard on saving the config file.
     1242 * \param *ConfigFileName name of file
     1243 * \param *configuration pointer to configuration structure with all the values
     1244 * \param *periode pointer to periodentafel structure with all the elements
     1245 * \param *molecules list of molecules structure with all the atoms and coordinates
     1246 */
     1247static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
     1248{
     1249  char filename[MAXSTRINGSIZE];
     1250  ofstream output;
     1251  molecule *mol = new molecule(periode);
     1252  mol->SetNameFromFilename(ConfigFileName);
     1253
     1254  if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
     1255    eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
     1256  }
     1257
     1258
     1259  // first save as PDB data
     1260  if (ConfigFileName != NULL)
     1261    strcpy(filename, ConfigFileName);
     1262  if (output == NULL)
     1263    strcpy(filename,"main_pcp_linux");
     1264  Log() << Verbose(0) << "Saving as pdb input ";
     1265  if (configuration->SavePDB(filename, molecules))
     1266    Log() << Verbose(0) << "done." << endl;
     1267  else
     1268    Log() << Verbose(0) << "failed." << endl;
     1269
     1270  // then save as tremolo data file
     1271  if (ConfigFileName != NULL)
     1272    strcpy(filename, ConfigFileName);
     1273  if (output == NULL)
     1274    strcpy(filename,"main_pcp_linux");
     1275  Log() << Verbose(0) << "Saving as tremolo data input ";
     1276  if (configuration->SaveTREMOLO(filename, molecules))
     1277    Log() << Verbose(0) << "done." << endl;
     1278  else
     1279    Log() << Verbose(0) << "failed." << endl;
     1280
     1281  // translate each to its center and merge all molecules in MoleculeListClass into this molecule
     1282  int N = molecules->ListOfMolecules.size();
     1283  int *src = new int[N];
     1284  N=0;
     1285  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     1286    src[N++] = (*ListRunner)->IndexNr;
     1287    (*ListRunner)->Translate(&(*ListRunner)->Center);
     1288  }
     1289  molecules->SimpleMultiAdd(mol, src, N);
     1290  delete[](src);
     1291
     1292  // ... and translate back
     1293  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     1294    (*ListRunner)->Center.Scale(-1.);
     1295    (*ListRunner)->Translate(&(*ListRunner)->Center);
     1296    (*ListRunner)->Center.Scale(-1.);
     1297  }
     1298
     1299  Log() << Verbose(0) << "Storing configuration ... " << endl;
     1300  // get correct valence orbitals
     1301  mol->CalculateOrbitals(*configuration);
     1302  configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
     1303  if (ConfigFileName != NULL) { // test the file name
     1304    strcpy(filename, ConfigFileName);
     1305    output.open(filename, ios::trunc);
     1306  } else if (strlen(configuration->configname) != 0) {
     1307    strcpy(filename, configuration->configname);
     1308    output.open(configuration->configname, ios::trunc);
     1309    } else {
     1310      strcpy(filename, DEFAULTCONFIG);
     1311      output.open(DEFAULTCONFIG, ios::trunc);
     1312    }
     1313  output.close();
     1314  output.clear();
     1315  Log() << Verbose(0) << "Saving of config file ";
     1316  if (configuration->Save(filename, periode, mol))
     1317    Log() << Verbose(0) << "successful." << endl;
     1318  else
     1319    Log() << Verbose(0) << "failed." << endl;
     1320
     1321  // and save to xyz file
     1322  if (ConfigFileName != NULL) {
     1323    strcpy(filename, ConfigFileName);
     1324    strcat(filename, ".xyz");
     1325    output.open(filename, ios::trunc);
     1326  }
     1327  if (output == NULL) {
     1328    strcpy(filename,"main_pcp_linux");
     1329    strcat(filename, ".xyz");
     1330    output.open(filename, ios::trunc);
     1331  }
     1332  Log() << Verbose(0) << "Saving of XYZ file ";
     1333  if (mol->MDSteps <= 1) {
     1334    if (mol->OutputXYZ(&output))
     1335      Log() << Verbose(0) << "successful." << endl;
     1336    else
     1337      Log() << Verbose(0) << "failed." << endl;
     1338  } else {
     1339    if (mol->OutputTrajectoriesXYZ(&output))
     1340      Log() << Verbose(0) << "successful." << endl;
     1341    else
     1342      Log() << Verbose(0) << "failed." << endl;
     1343  }
     1344  output.close();
     1345  output.clear();
     1346
     1347  // and save as MPQC configuration
     1348  if (ConfigFileName != NULL)
     1349    strcpy(filename, ConfigFileName);
     1350  if (output == NULL)
     1351    strcpy(filename,"main_pcp_linux");
     1352  Log() << Verbose(0) << "Saving as mpqc input ";
     1353  if (configuration->SaveMPQC(filename, mol))
     1354    Log() << Verbose(0) << "done." << endl;
     1355  else
     1356    Log() << Verbose(0) << "failed." << endl;
     1357
     1358  if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
     1359    eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
     1360  }
     1361
     1362  delete(mol);
     1363};
    761364
    771365/** Parses the command line options.
    … …  
    851373 * \return exit code (0 - successful, all else - something's wrong)
    861374 */
    87 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,\
    88                                    config& configuration, char *&ConfigFileName)
     1375static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName)
    891376{
    901377  Vector x,y,z,n;  // coordinates for absolute point in cell volume
    … …  
    1021389  int argptr;
    1031390  molecule *mol = NULL;
    104   string BondGraphFileName("");
     1391  string BondGraphFileName("\n");
    1051392  int verbosity = 0;
    1061393  strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);
    … …  
    1251412            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;
    1261413            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    127             Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
     1414            Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
    1281415            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    1291416            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    1301417            Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    1311418            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
    132             Log() << Verbose(0) << "\t-f/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;
     1419            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;
     1420            Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
    1331421            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
    1341422            Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
     1423            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
    1351424            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;
    1361425            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    2561545       mol = new molecule(periode);
    2571546       mol->ActiveFlag = true;
     1547       if (ConfigFileName != NULL)
     1548         mol->SetNameFromFilename(ConfigFileName);
    2581549       molecules->insert(mol);
     1550     }
     1551     if (configuration.BG == NULL) {
     1552       configuration.BG = new BondGraph(configuration.GetIsAngstroem());
     1553       if ((!BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) {
     1554         Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
     1555       } else {
     1556         eLog() << Verbose(1) << "Bond length table loading failed." << endl;
     1557       }
    2591558     }
    2601559
    … …  
    3581657              //argptr+=1;
    3591658              break;
     1659            case 'I':
     1660              Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;
     1661              // @TODO rather do the dissection afterwards
     1662              molecules->DissectMoleculeIntoConnectedSubgraphs(periode, &configuration);
     1663              mol = NULL;
     1664              if (molecules->ListOfMolecules.size() != 0) {
     1665                for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1666                  if ((*ListRunner)->ActiveFlag) {
     1667                    mol = *ListRunner;
     1668                    break;
     1669                  }
     1670              }
     1671              if (mol == NULL) {
     1672                mol = *(molecules->ListOfMolecules.begin());
     1673                mol->ActiveFlag = true;
     1674              }
     1675              break;
    3601676            case 'C':
    3611677              if (ExitFlag == 0) ExitFlag = 1;
    … …  
    3651681                performCriticalExit();
    3661682              } else {
    367                 SaveFlag = false;
    3681683                ofstream output(argv[argptr+1]);
    3691684                ofstream binoutput(argv[argptr+2]);
    … …  
    3851700                counter = 0;
    3861701                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    387                   Actives[counter] = (*BigFinder)->ActiveFlag;
     1702                  Actives[counter++] = (*BigFinder)->ActiveFlag;
    3881703                  (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
    3891704                }
    … …  
    3931708                int ranges[NDIM] = {1,1,1};
    3941709                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
     1710                //OutputCorrelationToSurface(&output, surfacemap);
    3951711                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
    3961712                OutputCorrelation ( &binoutput, binmap );
    … …  
    3981714                binoutput.close();
    3991715                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
    400                   (*BigFinder)->ActiveFlag = Actives[counter];
     1716                  (*BigFinder)->ActiveFlag = Actives[counter++];
    4011717                Free(&Actives);
    4021718                delete(LCList);
    … …  
    4211737            case 'F':
    4221738              if (ExitFlag == 0) ExitFlag = 1;
    423               if (argptr+5 >=argc) {
     1739              if (argptr+6 >=argc) {
    4241740                ExitFlag = 255;
    425                 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
     1741                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
    4261742                performCriticalExit();
    4271743              } else {
    … …  
    4291745                Log() << Verbose(1) << "Filling Box with water molecules." << endl;
    4301746                // construct water molecule
    431                 molecule *filler = new molecule(periode);;
     1747                molecule *filler = new molecule(periode);
    4321748                molecule *Filling = NULL;
    4331749                atom *second = NULL, *third = NULL;
     1750//                first = new atom();
     1751//                first->type = periode->FindElement(5);
     1752//                first->x.Zero();
     1753//                filler->AddAtom(first);
    4341754                first = new atom();
    4351755                first->type = periode->FindElement(1);
    … …  
    4501770                for (int i=0;i<NDIM;i++)
    4511771                  distance[i] = atof(argv[argptr+i]);
    452                 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
     1772                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
    4531773                if (Filling != NULL) {
     1774                  Filling->ActiveFlag = false;
    4541775                  molecules->insert(Filling);
    4551776                }
    … …  
    4981819                start = clock();
    4991820                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
    500                 FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]);
     1821                if (!FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]))
     1822                  ExitFlag = 255;
    5011823                //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
    5021824                end = clock();
    5031825                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    5041826                delete(LCList);
     1827                delete(T);
    5051828                argptr+=2;
    5061829              }
    … …  
    7982121                if (volume != -1)
    7992122                  ExitFlag = 255;
    800                   eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
     2123                  eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;
    8012124                  performCriticalExit();
    8022125              } else {
    … …  
    8872210    } while (argptr < argc);
    8882211    if (SaveFlag)
    889       configuration.SaveAll(ConfigFileName, periode, molecules);
     2212      SaveConfig(ConfigFileName, &configuration, periode, molecules);
    8902213  } else {  // no arguments, hence scan the elements db
    8912214    if (periode->LoadPeriodentafel(configuration.databasepath))
    … …  
    8982221};
    8992222
    900 /***************************************** Functions used to build all menus **********************/
    901 
    902 void populateEditMoleculesMenu(Menu* editMoleculesMenu,MoleculeListClass *molecules, config *configuration, periodentafel *periode){
    903   // build the EditMoleculesMenu
    904   Action *createMoleculeAction = new MethodAction("createMoleculeAction",boost::bind(&MoleculeListClass::createNewMolecule,molecules,periode));
    905   new ActionMenuItem('c',"create new molecule",editMoleculesMenu,createMoleculeAction);
    906 
    907   Action *loadMoleculeAction = new MethodAction("loadMoleculeAction",boost::bind(&MoleculeListClass::loadFromXYZ,molecules,periode));
    908   new ActionMenuItem('l',"load molecule from xyz file",editMoleculesMenu,loadMoleculeAction);
    909 
    910   Action *changeFilenameAction = new MethodAction("changeFilenameAction",boost::bind(&MoleculeListClass::changeName,molecules));
    911   new ActionMenuItem('n',"change molecule's name",editMoleculesMenu,changeFilenameAction);
    912 
    913   Action *giveFilenameAction = new MethodAction("giveFilenameAction",boost::bind(&MoleculeListClass::setMoleculeFilename,molecules));
    914   new ActionMenuItem('N',"give molecules filename",editMoleculesMenu,giveFilenameAction);
    915 
    916   Action *parseAtomsAction = new MethodAction("parseAtomsAction",boost::bind(&MoleculeListClass::parseXYZIntoMolecule,molecules));
    917   new ActionMenuItem('p',"parse atoms in xyz file into molecule",editMoleculesMenu,parseAtomsAction);
    918 
    919   Action *eraseMoleculeAction = new MethodAction("eraseMoleculeAction",boost::bind(&MoleculeListClass::eraseMolecule,molecules));
    920   new ActionMenuItem('r',"remove a molecule",editMoleculesMenu,eraseMoleculeAction);
    921 }
    922 
    923 
    9242223/********************************************** Main routine **************************************/
    9252224
    9262225int main(int argc, char **argv)
    9272226{
    928   periodentafel *periode = new periodentafel;
    929     MoleculeListClass *molecules = new MoleculeListClass;
    930     molecule *mol = NULL;
    931     config *configuration = new config;
    932     Vector x, y, z, n;
    933     ifstream test;
    934     ofstream output;
    935     string line;
    936     char *ConfigFileName = NULL;
    937     int j;
    938     setVerbosity(0);
    939     /* structure of ParseCommandLineOptions will be refactored later */
    940     j = ParseCommandLineOptions(argc, argv, molecules, periode, *configuration, ConfigFileName);
    941     switch (j){
    942         case 255:
    943         case 2:
    944         case 1:
    945             delete (molecules);
    946             delete (periode);
    947             delete (configuration);
    948             Log() << Verbose(0) << "Maximum of allocated memory: " << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
    949             Log() << Verbose(0) << "Remaining non-freed memory: " << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
    950             MemoryUsageObserver::getInstance()->purgeInstance();
    951             logger::purgeInstance();
    952             errorLogger::purgeInstance();
    953             return (j == 1 ? 0 : j);
    954         default:
    955             break;
     2227  periodentafel *periode = new periodentafel; // and a period table of all elements
     2228  MoleculeListClass *molecules = new MoleculeListClass;  // list of all molecules
     2229  molecule *mol = NULL;
     2230  config *configuration = new config;
     2231  char choice;  // menu choice char
     2232  Vector x,y,z,n;  // coordinates for absolute point in cell volume
     2233  ifstream test;
     2234  ofstream output;
     2235  string line;
     2236  char *ConfigFileName = NULL;
     2237  int j;
     2238
     2239  cout << ESPACKVersion << endl;
     2240
     2241  setVerbosity(0);
     2242
     2243  // =========================== PARSE COMMAND LINE OPTIONS ====================================
     2244  j = ParseCommandLineOptions(argc, argv, molecules, periode, *configuration, ConfigFileName);
     2245  switch(j) {
     2246    case 255:  // something went wrong
     2247    case 2:  // just for -f option
     2248    case 1:  // just for -v and -h options
     2249      delete(molecules); // also free's all molecules contained
     2250      delete(periode);
     2251      delete(configuration);
     2252      Log() << Verbose(0) <<  "Maximum of allocated memory: "
     2253        << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
     2254      Log() << Verbose(0) <<  "Remaining non-freed memory: "
     2255        << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
     2256      MemoryUsageObserver::getInstance()->purgeInstance();
     2257      logger::purgeInstance();
     2258      errorLogger::purgeInstance();
     2259     return (j == 1 ? 0 : j);
     2260    default:
     2261      break;
     2262  }
     2263
     2264  // General stuff
     2265  if (molecules->ListOfMolecules.size() == 0) {
     2266    mol = new molecule(periode);
     2267    if (mol->cell_size[0] == 0.) {
     2268      Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
     2269      for (int i=0;i<6;i++) {
     2270        Log() << Verbose(1) << "Cell size" << i << ": ";
     2271        cin >> mol->cell_size[i];
     2272      }
    9562273    }
    957     if(molecules->ListOfMolecules.size() == 0){
    958         mol = new molecule(periode);
    959         if(mol->cell_size[0] == 0.){
    960             Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
    961             for(int i = 0;i < 6;i++){
    962                 Log() << Verbose(1) << "Cell size" << i << ": ";
    963                 cin >> mol->cell_size[i];
    964             }
     2274    mol->ActiveFlag = true;
     2275    molecules->insert(mol);
     2276  }
     2277
     2278  // =========================== START INTERACTIVE SESSION ====================================
     2279
     2280  // now the main construction loop
     2281  Log() << Verbose(0) << endl << "Now comes the real construction..." << endl;
     2282  do {
     2283    Log() << Verbose(0) << endl << endl;
     2284    Log() << Verbose(0) << "============Molecule list=======================" << endl;
     2285    molecules->Enumerate((ofstream *)&cout);
     2286    Log() << Verbose(0) << "============Menu===============================" << endl;
     2287    Log() << Verbose(0) << "a - set molecule (in)active" << endl;
     2288    Log() << Verbose(0) << "e - edit molecules (load, parse, save)" << endl;
     2289    Log() << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;
     2290    Log() << Verbose(0) << "M - Merge molecules" << endl;
     2291    Log() << Verbose(0) << "m - manipulate atoms" << endl;
     2292    Log() << Verbose(0) << "-----------------------------------------------" << endl;
     2293    Log() << Verbose(0) << "c - edit the current configuration" << endl;
     2294    Log() << Verbose(0) << "-----------------------------------------------" << endl;
     2295    Log() << Verbose(0) << "s - save current setup to config file" << endl;
     2296    Log() << Verbose(0) << "T - call the current test routine" << endl;
     2297    Log() << Verbose(0) << "q - quit" << endl;
     2298    Log() << Verbose(0) << "===============================================" << endl;
     2299    Log() << Verbose(0) << "Input: ";
     2300    cin >> choice;
     2301
     2302    switch (choice) {
     2303      case 'a':  // (in)activate molecule
     2304        {
     2305          Log() << Verbose(0) << "Enter index of molecule: ";
     2306          cin >> j;
     2307          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     2308            if ((*ListRunner)->IndexNr == j)
     2309              (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
    9652310        }
    966 
    967         mol->ActiveFlag = true;
    968         molecules->insert(mol);
    969     }
    970 
    971     {
    972       menuPopulaters populaters;
    973       populaters.MakeEditMoleculesMenu = populateEditMoleculesMenu;
    974 
    975       UIFactory::makeUserInterface(UIFactory::Text);
    976       MainWindow *mainWindow = UIFactory::get()->makeMainWindow(populaters,molecules, configuration, periode, ConfigFileName);
    977       mainWindow->display();
    978       delete mainWindow;
    979     }
    980 
    981     if(periode->StorePeriodentafel(configuration->databasepath))
    982         Log() << Verbose(0) << "Saving of elements.db successful." << endl;
    983 
    984     else
    985         Log() << Verbose(0) << "Saving of elements.db failed." << endl;
    986 
    987     delete (molecules);
    988     delete(periode);
     2311        break;
     2312
     2313      case 'c': // edit each field of the configuration
     2314       configuration->Edit();
     2315       break;
     2316
     2317      case 'e': // create molecule
     2318        EditMolecules(periode, molecules);
     2319        break;
     2320
     2321      case 'g': // manipulate molecules
     2322        ManipulateMolecules(periode, molecules, configuration);
     2323        break;
     2324
     2325      case 'M':  // merge molecules
     2326        MergeMolecules(periode, molecules);
     2327        break;
     2328
     2329      case 'm': // manipulate atoms
     2330        ManipulateAtoms(periode, molecules, configuration);
     2331        break;
     2332
     2333      case 'q': // quit
     2334        break;
     2335
     2336      case 's': // save to config file
     2337        SaveConfig(ConfigFileName, configuration, periode, molecules);
     2338        break;
     2339
     2340      case 'T':
     2341        testroutine(molecules);
     2342        break;
     2343
     2344      default:
     2345        break;
     2346    };
     2347  } while (choice != 'q');
     2348
     2349  // save element data base
     2350  if (periode->StorePeriodentafel(configuration->databasepath)) //ElementsFileName
     2351    Log() << Verbose(0) << "Saving of elements.db successful." << endl;
     2352  else
     2353    Log() << Verbose(0) << "Saving of elements.db failed." << endl;
     2354
     2355  delete(molecules); // also free's all molecules contained
     2356  delete(periode);
    9892357  delete(configuration);
    990 
    991 
    9922358
    9932359  Log() << Verbose(0) <<  "Maximum of allocated memory: "
    … …  
    9982364  logger::purgeInstance();
    9992365  errorLogger::purgeInstance();
    1000   UIFactory::purgeInstance();
    1001   ActionRegistry::purgeRegistry();
     2366
    10022367  return (0);
    10032368}
Note: See TracChangeset for help on using the changeset viewer.