Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    r244a84 ra6f180  
    4848
    4949
     50#include <boost/bind.hpp>
     51
    5052using namespace std;
    51 
    52 #include <cstring>
    5353
    5454#include "analysis_correlation.hpp"
    … …  
    6767#include "molecule.hpp"
    6868#include "periodentafel.hpp"
    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  */
    77 static 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  */
    275 static 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  */
    333 static 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  */
    401 static 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  */
    448 static 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  */
    512 static 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  */
    632 static 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  */
    655 static 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  */
    774 static 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  */
    933 static 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  */
    1049 static 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  */
    1163 static 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  */
    1247 static 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 };
     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#include "Actions/small_actions.hpp"
     76
    136477
    136578/** Parses the command line options.
    … …  
    137386 * \return exit code (0 - successful, all else - something's wrong)
    137487 */
    1375 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName)
     88static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,\
     89                                   config& configuration, char *&ConfigFileName)
    137690{
    137791  Vector x,y,z,n;  // coordinates for absolute point in cell volume
    … …  
    1389103  int argptr;
    1390104  molecule *mol = NULL;
    1391   string BondGraphFileName("\n");
     105  string BondGraphFileName("");
    1392106  int verbosity = 0;
    1393107  strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);
    … …  
    1412126            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;
    1413127            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    1414             Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
     128            Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
    1415129            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    1416130            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    1417131            Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    1418132            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << 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;
     133            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;
    1421134            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
    1422135            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;
    1424136            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;
    1425137            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    1545257       mol = new molecule(periode);
    1546258       mol->ActiveFlag = true;
    1547        if (ConfigFileName != NULL)
    1548          mol->SetNameFromFilename(ConfigFileName);
    1549259       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        }
    1558260     }
    1559261
    … …  
    1657359              //argptr+=1;
    1658360              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;
    1676361            case 'C':
    1677362              if (ExitFlag == 0) ExitFlag = 1;
    … …  
    1681366                performCriticalExit();
    1682367              } else {
     368                SaveFlag = false;
    1683369                ofstream output(argv[argptr+1]);
    1684370                ofstream binoutput(argv[argptr+2]);
    … …  
    1700386                counter = 0;
    1701387                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1702                   Actives[counter++] = (*BigFinder)->ActiveFlag;
     388                  Actives[counter] = (*BigFinder)->ActiveFlag;
    1703389                  (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
    1704390                }
    … …  
    1708394                int ranges[NDIM] = {1,1,1};
    1709395                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
    1710                 //OutputCorrelationToSurface(&output, surfacemap);
    1711396                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
    1712397                OutputCorrelation ( &binoutput, binmap );
    … …  
    1714399                binoutput.close();
    1715400                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
    1716                   (*BigFinder)->ActiveFlag = Actives[counter++];
     401                  (*BigFinder)->ActiveFlag = Actives[counter];
    1717402                Free(&Actives);
    1718403                delete(LCList);
    … …  
    1737422            case 'F':
    1738423              if (ExitFlag == 0) ExitFlag = 1;
    1739               if (argptr+6 >=argc) {
    1740                 ExitFlag = 255;
    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;
     424              if (argptr+5 >=argc) {
     425                ExitFlag = 255;
     426                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;
    1742427                performCriticalExit();
    1743428              } else {
    … …  
    1745430                Log() << Verbose(1) << "Filling Box with water molecules." << endl;
    1746431                // construct water molecule
    1747                 molecule *filler = new molecule(periode);
     432                molecule *filler = new molecule(periode);;
    1748433                molecule *Filling = NULL;
    1749434                atom *second = NULL, *third = NULL;
    1750 //                first = new atom();
    1751 //                first->type = periode->FindElement(5);
    1752 //                first->x.Zero();
    1753 //                filler->AddAtom(first);
    1754435                first = new atom();
    1755436                first->type = periode->FindElement(1);
    … …  
    1770451                for (int i=0;i<NDIM;i++)
    1771452                  distance[i] = atof(argv[argptr+i]);
    1772                 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
     453                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
    1773454                if (Filling != NULL) {
    1774                   Filling->ActiveFlag = false;
    1775455                  molecules->insert(Filling);
    1776456                }
    … …  
    1819499                start = clock();
    1820500                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
    1821                 if (!FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]))
    1822                   ExitFlag = 255;
     501                FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]);
    1823502                //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
    1824503                end = clock();
    1825504                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    1826505                delete(LCList);
    1827                 delete(T);
    1828506                argptr+=2;
    1829507              }
    … …  
    2121799                if (volume != -1)
    2122800                  ExitFlag = 255;
    2123                   eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;
     801                  eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
    2124802                  performCriticalExit();
    2125803              } else {
    … …  
    2210888    } while (argptr < argc);
    2211889    if (SaveFlag)
    2212       SaveConfig(ConfigFileName, &configuration, periode, molecules);
     890      configuration.SaveAll(ConfigFileName, periode, molecules);
    2213891  } else {  // no arguments, hence scan the elements db
    2214892    if (periode->LoadPeriodentafel(configuration.databasepath))
    … …  
    2221899};
    2222900
     901/***************************************** Functions used to build all menus **********************/
     902
     903void populateEditMoleculesMenu(Menu* editMoleculesMenu,MoleculeListClass *molecules, config *configuration, periodentafel *periode){
     904  // build the EditMoleculesMenu
     905  Action *createMoleculeAction = new MethodAction("createMoleculeAction",boost::bind(&MoleculeListClass::createNewMolecule,molecules,periode));
     906  new ActionMenuItem('c',"create new molecule",editMoleculesMenu,createMoleculeAction);
     907
     908  Action *loadMoleculeAction = new MethodAction("loadMoleculeAction",boost::bind(&MoleculeListClass::loadFromXYZ,molecules,periode));
     909  new ActionMenuItem('l',"load molecule from xyz file",editMoleculesMenu,loadMoleculeAction);
     910
     911  Action *changeFilenameAction = new ChangeMoleculeNameAction(molecules);
     912  new ActionMenuItem('n',"change molecule's name",editMoleculesMenu,changeFilenameAction);
     913
     914  Action *giveFilenameAction = new MethodAction("giveFilenameAction",boost::bind(&MoleculeListClass::setMoleculeFilename,molecules));
     915  new ActionMenuItem('N',"give molecules filename",editMoleculesMenu,giveFilenameAction);
     916
     917  Action *parseAtomsAction = new MethodAction("parseAtomsAction",boost::bind(&MoleculeListClass::parseXYZIntoMolecule,molecules));
     918  new ActionMenuItem('p',"parse atoms in xyz file into molecule",editMoleculesMenu,parseAtomsAction);
     919
     920  Action *eraseMoleculeAction = new MethodAction("eraseMoleculeAction",boost::bind(&MoleculeListClass::eraseMolecule,molecules));
     921  new ActionMenuItem('r',"remove a molecule",editMoleculesMenu,eraseMoleculeAction);
     922}
     923
     924
    2223925/********************************************** Main routine **************************************/
    2224926
    2225927int main(int argc, char **argv)
    2226928{
    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       }
     929  periodentafel *periode = new periodentafel;
     930    MoleculeListClass *molecules = new MoleculeListClass;
     931    molecule *mol = NULL;
     932    config *configuration = new config;
     933    Vector x, y, z, n;
     934    ifstream test;
     935    ofstream output;
     936    string line;
     937    char *ConfigFileName = NULL;
     938    int j;
     939    setVerbosity(0);
     940    /* structure of ParseCommandLineOptions will be refactored later */
     941    j = ParseCommandLineOptions(argc, argv, molecules, periode, *configuration, ConfigFileName);
     942    switch (j){
     943        case 255:
     944        case 2:
     945        case 1:
     946            delete (molecules);
     947            delete (periode);
     948            delete (configuration);
     949            Log() << Verbose(0) << "Maximum of allocated memory: " << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
     950            Log() << Verbose(0) << "Remaining non-freed memory: " << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
     951            MemoryUsageObserver::getInstance()->purgeInstance();
     952            logger::purgeInstance();
     953            errorLogger::purgeInstance();
     954            return (j == 1 ? 0 : j);
     955        default:
     956            break;
    2273957    }
    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;
     958    if(molecules->ListOfMolecules.size() == 0){
     959        mol = new molecule(periode);
     960        if(mol->cell_size[0] == 0.){
     961            Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
     962            for(int i = 0;i < 6;i++){
     963                Log() << Verbose(1) << "Cell size" << i << ": ";
     964                cin >> mol->cell_size[i];
     965            }
    2310966        }
    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);
     967
     968        mol->ActiveFlag = true;
     969        molecules->insert(mol);
     970    }
     971
     972    {
     973      menuPopulaters populaters;
     974      populaters.MakeEditMoleculesMenu = populateEditMoleculesMenu;
     975
     976#ifdef USE_GUI_QT
     977      UIFactory::makeUserInterface(UIFactory::QT4);
     978#else
     979      UIFactory::makeUserInterface(UIFactory::Text);
     980#endif
     981      MainWindow *mainWindow = UIFactory::get()->makeMainWindow(populaters,molecules, configuration, periode, ConfigFileName);
     982      mainWindow->display();
     983      delete mainWindow;
     984    }
     985
     986    if(periode->StorePeriodentafel(configuration->databasepath))
     987        Log() << Verbose(0) << "Saving of elements.db successful." << endl;
     988
     989    else
     990        Log() << Verbose(0) << "Saving of elements.db failed." << endl;
     991
     992    delete (molecules);
     993    delete(periode);
    2357994  delete(configuration);
     995
     996
    2358997
    2359998  Log() << Verbose(0) <<  "Maximum of allocated memory: "
    … …  
    23641003  logger::purgeInstance();
    23651004  errorLogger::purgeInstance();
    2366 
     1005  UIFactory::purgeInstance();
     1006  ActionRegistry::purgeRegistry();
    23671007  return (0);
    23681008}
Note: See TracChangeset for help on using the changeset viewer.