Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

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