Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    rda3024e ra356f2  
    11/** \file builder.cpp
    22 *
    3  * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
    4  * The output is the complete configuration file for PCP for direct use.
    5  * Features:
    6  * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
    7  * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
     3 *  date: Jan 1, 2007
     4 *  author: heber
    85 *
    96 */
    107
    11 /*! \mainpage Molecuilder - a molecular set builder
     8/*! \mainpage MoleCuilder - a molecular set builder
    129 *
    13  * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
     10 * This introductory shall briefly make acquainted with the program, helping in installing and a first run.
    1411 *
    1512 * \section about About the Program
    1613 *
    17  *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
    18  *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
    19  *  already constructed atoms.
     14 *  MoleCuilder is a program, written entirely in C++, that enables the construction of a coordinate set for the
     15 *  atoms making up an molecule. It allows for both building of simple molecules by adding atom-wise giving bond
     16 *  angles and distances or absolute coordinates, but also using them as templates. Regions can be specified and
     17 *  ordered to be filled with a molecule in a certain manner. Greater conglomerations of molecules can be tesselated
     18 *  and recognized as a region themselves to be subsequently surrounded by other (surface solvated) molecules.
     19 *  In the end, MoleCuilder allows the construction of arbitrary nano structures, whether they be crystalline or
     20 *  amorphic in nature.
    2021 *
    21  *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
    22  *  molecular dynamics implementation.
    2322 *
    2423 * \section install Installation
    … …  
    3231 *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
    3332 *  -# make doxygen-doc: Creates these html pages out of the documented source
     33 *  -# make check: Runs an extensive set of unit tests and a testsuite which also gives a good overview on the set of
     34 *                 functions.
    3435 *
    3536 * \section run Running
    … …  
    3738 *  The program can be executed by running: ./molecuilder
    3839 *
    39  *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
    40  *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
    41  *  later re-execution.
     40 *  MoleCuilder has three interfaces at your disposal:
     41 *  -# Textmenu: A simple interactive console-based menu, where awaits your choices and inputs in order to set atoms
     42 *               as you like
     43 *  -# CommandLineUI: Every command can also be chained up as a sequence of actions on the command line to be executed
     44 *               with any user interaction.
     45 *  -# GraphicalUI: A graphical user interface that also display the molecular structure being built and lots of other
     46 *               informations to ease the construction of bigger geometries.
    4247 *
    43  * \section ref References
    44  *
    45  *  For the special configuration file format, see the documentation of pcp.
     48 *  The supported output formats right now are:
     49 *  -# mpqc: Configuration files of the Massively Parallel Quantum Chemistry package (Sandia labs)
     50 *  -# pcp: Configuration files of the Parallel Car-Parrinello program (Institute for Numerical Simulation)
     51 *  -# tremolo: Configuration files of TREMOLO (Institute for Numerical Simulation)
     52 *  -# xyz: the most basic format for the 3d arrangement of atoms consisting of a list of element and 3 coordinates.
    4653 *
    4754 */
    … …  
    4956#include "Helpers/MemDebug.hpp"
    5057
    51 #include <boost/bind.hpp>
    52 
    53 using namespace std;
    54 
    55 #include <cstring>
    56 #include <cstdlib>
    57 
    58 #include "analysis_bonds.hpp"
    59 #include "analysis_correlation.hpp"
    60 #include "atom.hpp"
    61 #include "bond.hpp"
    6258#include "bondgraph.hpp"
    63 #include "boundary.hpp"
    6459#include "CommandLineParser.hpp"
    6560#include "config.hpp"
    66 #include "element.hpp"
    67 #include "ellipsoid.hpp"
    68 #include "helpers.hpp"
    69 #include "leastsquaremin.hpp"
    70 #include "linkedcell.hpp"
    7161#include "log.hpp"
    72 #include "memoryusageobserver.hpp"
    73 #include "molecule.hpp"
    74 #include "periodentafel.hpp"
     62#include "verbose.hpp"
     63#include "World.hpp"
     64
     65#include "Actions/ActionRegistry.hpp"
     66#include "Actions/ActionHistory.hpp"
     67#include "Actions/MapOfActions.hpp"
     68
     69#include "Parser/ChangeTracker.hpp"
     70#include "Parser/FormatParserStorage.hpp"
     71
    7572#include "UIElements/UIFactory.hpp"
    7673#include "UIElements/TextUI/TextUIFactory.hpp"
    … …  
    7875#include "UIElements/MainWindow.hpp"
    7976#include "UIElements/Dialog.hpp"
    80 #include "Menu/ActionMenuItem.hpp"
    81 #include "Actions/ActionRegistry.hpp"
    82 #include "Actions/ActionHistory.hpp"
    83 #include "Actions/MapOfActions.hpp"
    84 #include "Actions/MethodAction.hpp"
    85 #include "Actions/MoleculeAction/ChangeNameAction.hpp"
    86 #include "World.hpp"
     77
    8778#include "version.h"
    88 #include "World.hpp"
    8979
    90 
    91 /********************************************* Subsubmenu routine ************************************/
    92 #if 0
    93 /** Submenu for adding atoms to the molecule.
    94  * \param *periode periodentafel
    95  * \param *molecule molecules with atoms
    96  */
    97 static void AddAtoms(periodentafel *periode, molecule *mol)
    98 {
    99   atom *first, *second, *third, *fourth;
    100   Vector **atoms;
    101   Vector x,y,z,n;  // coordinates for absolute point in cell volume
    102   double a,b,c;
    103   char choice;  // menu choice char
    104   bool valid;
    105 
    106   cout << Verbose(0) << "===========ADD ATOM============================" << endl;
    107   cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
    108   cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
    109   cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
    110   cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
    111   cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
    112   cout << Verbose(0) << "all else - go back" << endl;
    113   cout << Verbose(0) << "===============================================" << endl;
    114   cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
    115   cout << Verbose(0) << "INPUT: ";
    116   cin >> choice;
    117 
    118   switch (choice) {
    119     default:
    120       DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl);
    121       break;
    122       case 'a': // absolute coordinates of atom
    123         cout << Verbose(0) << "Enter absolute coordinates." << endl;
    124         first = new atom;
    125         first->x.AskPosition(World::getInstance().getDomain(), false);
    126         first->type = periode->AskElement();  // give type
    127         mol->AddAtom(first);  // add to molecule
    128         break;
    129 
    130       case 'b': // relative coordinates of atom wrt to reference point
    131         first = new atom;
    132         valid = true;
    133         do {
    134           if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
    135           cout << Verbose(0) << "Enter reference coordinates." << endl;
    136           x.AskPosition(World::getInstance().getDomain(), true);
    137           cout << Verbose(0) << "Enter relative coordinates." << endl;
    138           first->x.AskPosition(World::getInstance().getDomain(), false);
    139           first->x.AddVector((const Vector *)&x);
    140           cout << Verbose(0) << "\n";
    141         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    142         first->type = periode->AskElement();  // give type
    143         mol->AddAtom(first);  // add to molecule
    144         break;
    145 
    146       case 'c': // relative coordinates of atom wrt to already placed atom
    147         first = new atom;
    148         valid = true;
    149         do {
    150           if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
    151           second = mol->AskAtom("Enter atom number: ");
    152           DoLog(0) && (Log() << Verbose(0) << "Enter relative coordinates." << endl);
    153           first->x.AskPosition(World::getInstance().getDomain(), false);
    154           for (int i=NDIM;i--;) {
    155             first->x.x[i] += second->x.x[i];
    156           }
    157         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    158         first->type = periode->AskElement();  // give type
    159         mol->AddAtom(first);  // add to molecule
    160         break;
    161 
    162     case 'd': // two atoms, two angles and a distance
    163         first = new atom;
    164         valid = true;
    165         do {
    166           if (!valid) {
    167             DoeLog(2) && (eLog()<< Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl);
    168           }
    169           cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
    170           second = mol->AskAtom("Enter central atom: ");
    171           third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
    172           fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
    173           a = ask_value("Enter distance between central (first) and new atom: ");
    174           b = ask_value("Enter angle between new, first and second atom (degrees): ");
    175           b *= M_PI/180.;
    176           bound(&b, 0., 2.*M_PI);
    177           c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
    178           c *= M_PI/180.;
    179           bound(&c, -M_PI, M_PI);
    180           cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
    181 /*
    182           second->Output(1,1,(ofstream *)&cout);
    183           third->Output(1,2,(ofstream *)&cout);
    184           fourth->Output(1,3,(ofstream *)&cout);
    185           n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
    186           x.Copyvector(&second->x);
    187           x.SubtractVector(&third->x);
    188           x.Copyvector(&fourth->x);
    189           x.SubtractVector(&third->x);
    190 
    191           if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
    192          coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
    193             continue;
    194           }
    195           DoLog(0) && (Log() << Verbose(0) << "resulting relative coordinates: ");
    196           z.Output();
    197           DoLog(0) && (Log() << Verbose(0) << endl);
    198           */
    199           // calc axis vector
    200           x.CopyVector(&second->x);
    201           x.SubtractVector(&third->x);
    202           x.Normalize();
    203           Log() << Verbose(0) << "x: ",
    204           x.Output();
    205           DoLog(0) && (Log() << Verbose(0) << endl);
    206           z.MakeNormalVector(&second->x,&third->x,&fourth->x);
    207           Log() << Verbose(0) << "z: ",
    208           z.Output();
    209           DoLog(0) && (Log() << Verbose(0) << endl);
    210           y.MakeNormalVector(&x,&z);
    211           Log() << Verbose(0) << "y: ",
    212           y.Output();
    213           DoLog(0) && (Log() << Verbose(0) << endl);
    214 
    215           // rotate vector around first angle
    216           first->x.CopyVector(&x);
    217           first->x.RotateVector(&z,b - M_PI);
    218           Log() << Verbose(0) << "Rotated vector: ",
    219           first->x.Output();
    220           DoLog(0) && (Log() << Verbose(0) << endl);
    221           // remove the projection onto the rotation plane of the second angle
    222           n.CopyVector(&y);
    223           n.Scale(first->x.ScalarProduct(&y));
    224           Log() << Verbose(0) << "N1: ",
    225           n.Output();
    226           DoLog(0) && (Log() << Verbose(0) << endl);
    227           first->x.SubtractVector(&n);
    228           Log() << Verbose(0) << "Subtracted vector: ",
    229           first->x.Output();
    230           DoLog(0) && (Log() << Verbose(0) << endl);
    231           n.CopyVector(&z);
    232           n.Scale(first->x.ScalarProduct(&z));
    233           Log() << Verbose(0) << "N2: ",
    234           n.Output();
    235           DoLog(0) && (Log() << Verbose(0) << endl);
    236           first->x.SubtractVector(&n);
    237           Log() << Verbose(0) << "2nd subtracted vector: ",
    238           first->x.Output();
    239           DoLog(0) && (Log() << Verbose(0) << endl);
    240 
    241           // rotate another vector around second angle
    242           n.CopyVector(&y);
    243           n.RotateVector(&x,c - M_PI);
    244           Log() << Verbose(0) << "2nd Rotated vector: ",
    245           n.Output();
    246           DoLog(0) && (Log() << Verbose(0) << endl);
    247 
    248           // add the two linear independent vectors
    249           first->x.AddVector(&n);
    250           first->x.Normalize();
    251           first->x.Scale(a);
    252           first->x.AddVector(&second->x);
    253 
    254           DoLog(0) && (Log() << Verbose(0) << "resulting coordinates: ");
    255           first->x.Output();
    256           DoLog(0) && (Log() << Verbose(0) << endl);
    257         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    258         first->type = periode->AskElement();  // give type
    259         mol->AddAtom(first);  // add to molecule
    260         break;
    261 
    262       case 'e': // least square distance position to a set of atoms
    263         first = new atom;
    264         atoms = new (Vector*[128]);
    265         valid = true;
    266         for(int i=128;i--;)
    267           atoms[i] = NULL;
    268         int i=0, j=0;
    269         cout << Verbose(0) << "Now we need at least three molecules.\n";
    270         do {
    271           cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
    272           cin >> j;
    273           if (j != -1) {
    274             second = mol->FindAtom(j);
    275             atoms[i++] = &(second->x);
    276           }
    277         } while ((j != -1) && (i<128));
    278         if (i >= 2) {
    279           first->x.LSQdistance((const Vector **)atoms, i);
    280           first->x.Output();
    281           first->type = periode->AskElement();  // give type
    282           mol->AddAtom(first);  // add to molecule
    283         } else {
    284           delete first;
    285           cout << Verbose(0) << "Please enter at least two vectors!\n";
    286         }
    287         break;
    288   };
    289 };
    290 
    291 /** Submenu for centering the atoms in the molecule.
    292  * \param *mol molecule with all the atoms
    293  */
    294 static void CenterAtoms(molecule *mol)
    295 {
    296   Vector x, y, helper;
    297   char choice;  // menu choice char
    298 
    299   cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
    300   cout << Verbose(0) << " a - on origin" << endl;
    301   cout << Verbose(0) << " b - on center of gravity" << endl;
    302   cout << Verbose(0) << " c - within box with additional boundary" << endl;
    303   cout << Verbose(0) << " d - within given simulation box" << endl;
    304   cout << Verbose(0) << "all else - go back" << endl;
    305   cout << Verbose(0) << "===============================================" << endl;
    306   cout << Verbose(0) << "INPUT: ";
    307   cin >> choice;
    308 
    309   switch (choice) {
    310     default:
    311       cout << Verbose(0) << "Not a valid choice." << endl;
    312       break;
    313     case 'a':
    314       cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
    315       mol->CenterOrigin();
    316       break;
    317     case 'b':
    318       cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
    319       mol->CenterPeriodic();
    320       break;
    321     case 'c':
    322       cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
    323       for (int i=0;i<NDIM;i++) {
    324         cout << Verbose(0) << "Enter axis " << i << " boundary: ";
    325         cin >> y.x[i];
    326       }
    327       mol->CenterEdge(&x);  // make every coordinate positive
    328       mol->Center.AddVector(&y); // translate by boundary
    329       helper.CopyVector(&y);
    330       helper.Scale(2.);
    331       helper.AddVector(&x);
    332       mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
    333       break;
    334     case 'd':
    335       cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
    336       for (int i=0;i<NDIM;i++) {
    337         cout << Verbose(0) << "Enter axis " << i << " boundary: ";
    338         cin >> x.x[i];
    339       }
    340       // update Box of atoms by boundary
    341       mol->SetBoxDimension(&x);
    342       // center
    343       mol->CenterInBox();
    344       break;
    345   }
    346 };
    347 
    348 /** Submenu for aligning the atoms in the molecule.
    349  * \param *periode periodentafel
    350  * \param *mol molecule with all the atoms
    351  */
    352 static void AlignAtoms(periodentafel *periode, molecule *mol)
    353 {
    354   atom *first, *second, *third;
    355   Vector x,n;
    356   char choice;  // menu choice char
    357 
    358   cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
    359   cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
    360   cout << Verbose(0) << " b - state alignment vector" << endl;
    361   cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
    362   cout << Verbose(0) << " d - align automatically by least square fit" << endl;
    363   cout << Verbose(0) << "all else - go back" << endl;
    364   cout << Verbose(0) << "===============================================" << endl;
    365   cout << Verbose(0) << "INPUT: ";
    366   cin >> choice;
    367 
    368   switch (choice) {
    369     default:
    370     case 'a': // three atoms defining mirror plane
    371       first = mol->AskAtom("Enter first atom: ");
    372       second = mol->AskAtom("Enter second atom: ");
    373       third = mol->AskAtom("Enter third atom: ");
    374 
    375       n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
    376       break;
    377     case 'b': // normal vector of mirror plane
    378       cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    379       n.AskPosition(World::getInstance().getDomain(),0);
    380       n.Normalize();
    381       break;
    382     case 'c': // three atoms defining mirror plane
    383       first = mol->AskAtom("Enter first atom: ");
    384       second = mol->AskAtom("Enter second atom: ");
    385 
    386       n.CopyVector((const Vector *)&first->x);
    387       n.SubtractVector((const Vector *)&second->x);
    388       n.Normalize();
    389       break;
    390     case 'd':
    391       char shorthand[4];
    392       Vector a;
    393       struct lsq_params param;
    394       do {
    395         fprintf(stdout, "Enter the element of atoms to be chosen: ");
    396         fscanf(stdin, "%3s", shorthand);
    397       } while ((param.type = periode->FindElement(shorthand)) == NULL);
    398       cout << Verbose(0) << "Element is " << param.type->name << endl;
    399       mol->GetAlignvector(&param);
    400       for (int i=NDIM;i--;) {
    401         x.x[i] = gsl_vector_get(param.x,i);
    402         n.x[i] = gsl_vector_get(param.x,i+NDIM);
    403       }
    404       gsl_vector_free(param.x);
    405       cout << Verbose(0) << "Offset vector: ";
    406       x.Output();
    407       DoLog(0) && (Log() << Verbose(0) << endl);
    408       n.Normalize();
    409       break;
    410   };
    411   DoLog(0) && (Log() << Verbose(0) << "Alignment vector: ");
    412   n.Output();
    413   DoLog(0) && (Log() << Verbose(0) << endl);
    414   mol->Align(&n);
    415 };
    416 
    417 /** Submenu for mirroring the atoms in the molecule.
    418  * \param *mol molecule with all the atoms
    419  */
    420 static void MirrorAtoms(molecule *mol)
    421 {
    422   atom *first, *second, *third;
    423   Vector n;
    424   char choice;  // menu choice char
    425 
    426   DoLog(0) && (Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl);
    427   DoLog(0) && (Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl);
    428   DoLog(0) && (Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl);
    429   DoLog(0) && (Log() << Verbose(0) << " c - state two atoms in normal direction" << endl);
    430   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    431   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    432   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    433   cin >> choice;
    434 
    435   switch (choice) {
    436     default:
    437     case 'a': // three atoms defining mirror plane
    438       first = mol->AskAtom("Enter first atom: ");
    439       second = mol->AskAtom("Enter second atom: ");
    440       third = mol->AskAtom("Enter third atom: ");
    441 
    442       n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
    443       break;
    444     case 'b': // normal vector of mirror plane
    445       DoLog(0) && (Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl);
    446       n.AskPosition(World::getInstance().getDomain(),0);
    447       n.Normalize();
    448       break;
    449     case 'c': // three atoms defining mirror plane
    450       first = mol->AskAtom("Enter first atom: ");
    451       second = mol->AskAtom("Enter second atom: ");
    452 
    453       n.CopyVector((const Vector *)&first->x);
    454       n.SubtractVector((const Vector *)&second->x);
    455       n.Normalize();
    456       break;
    457   };
    458   DoLog(0) && (Log() << Verbose(0) << "Normal vector: ");
    459   n.Output();
    460   DoLog(0) && (Log() << Verbose(0) << endl);
    461   mol->Mirror((const Vector *)&n);
    462 };
    463 
    464 /** Submenu for removing the atoms from the molecule.
    465  * \param *mol molecule with all the atoms
    466  */
    467 static void RemoveAtoms(molecule *mol)
    468 {
    469   atom *first, *second;
    470   int axis;
    471   double tmp1, tmp2;
    472   char choice;  // menu choice char
    473 
    474   DoLog(0) && (Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl);
    475   DoLog(0) && (Log() << Verbose(0) << " a - state atom for removal by number" << endl);
    476   DoLog(0) && (Log() << Verbose(0) << " b - keep only in radius around atom" << endl);
    477   DoLog(0) && (Log() << Verbose(0) << " c - remove this with one axis greater value" << endl);
    478   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    479   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    480   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    481   cin >> choice;
    482 
    483   switch (choice) {
    484     default:
    485     case 'a':
    486       if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
    487         DoLog(1) && (Log() << Verbose(1) << "Atom removed." << endl);
    488       else
    489         DoLog(1) && (Log() << Verbose(1) << "Atom not found." << endl);
    490       break;
    491     case 'b':
    492       second = mol->AskAtom("Enter number of atom as reference point: ");
    493       DoLog(0) && (Log() << Verbose(0) << "Enter radius: ");
    494       cin >> tmp1;
    495       first = mol->start;
    496       second = first->next;
    497       while(second != mol->end) {
    498         first = second;
    499         second = first->next;
    500         if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
    501           mol->RemoveAtom(first);
    502       }
    503       break;
    504     case 'c':
    505       DoLog(0) && (Log() << Verbose(0) << "Which axis is it: ");
    506       cin >> axis;
    507       DoLog(0) && (Log() << Verbose(0) << "Lower boundary: ");
    508       cin >> tmp1;
    509       DoLog(0) && (Log() << Verbose(0) << "Upper boundary: ");
    510       cin >> tmp2;
    511       first = mol->start;
    512       second = first->next;
    513       while(second != mol->end) {
    514         first = second;
    515         second = first->next;
    516         if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
    517           //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
    518           mol->RemoveAtom(first);
    519         }
    520       }
    521       break;
    522   };
    523   //mol->Output();
    524   choice = 'r';
    525 };
    526 
    527 /** Submenu for measuring out the atoms in the molecule.
    528  * \param *periode periodentafel
    529  * \param *mol molecule with all the atoms
    530  */
    531 static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
    532 {
    533   atom *first, *second, *third;
    534   Vector x,y;
    535   double min[256], tmp1, tmp2, tmp3;
    536   int Z;
    537   char choice;  // menu choice char
    538 
    539   DoLog(0) && (Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl);
    540   DoLog(0) && (Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl);
    541   DoLog(0) && (Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl);
    542   DoLog(0) && (Log() << Verbose(0) << " c - calculate bond angle" << endl);
    543   DoLog(0) && (Log() << Verbose(0) << " d - calculate principal axis of the system" << endl);
    544   DoLog(0) && (Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl);
    545   DoLog(0) && (Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl);
    546   DoLog(0) && (Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl);
    547   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    548   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    549   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    550   cin >> choice;
    551 
    552   switch(choice) {
    553     default:
    554       DoLog(1) && (Log() << Verbose(1) << "Not a valid choice." << endl);
    555       break;
    556     case 'a':
    557       first = mol->AskAtom("Enter first atom: ");
    558       for (int i=MAX_ELEMENTS;i--;)
    559         min[i] = 0.;
    560 
    561       second = mol->start;
    562       while ((second->next != mol->end)) {
    563         second = second->next; // advance
    564         Z = second->type->Z;
    565         tmp1 = 0.;
    566         if (first != second) {
    567           x.CopyVector((const Vector *)&first->x);
    568           x.SubtractVector((const Vector *)&second->x);
    569           tmp1 = x.Norm();
    570         }
    571         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
    572         //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
    573       }
    574       for (int i=MAX_ELEMENTS;i--;)
    575         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;
    576       break;
    577 
    578     case 'b':
    579       first = mol->AskAtom("Enter first atom: ");
    580       second = mol->AskAtom("Enter second atom: ");
    581       for (int i=NDIM;i--;)
    582         min[i] = 0.;
    583       x.CopyVector((const Vector *)&first->x);
    584       x.SubtractVector((const Vector *)&second->x);
    585       tmp1 = x.Norm();
    586       DoLog(1) && (Log() << Verbose(1) << "Distance vector is ");
    587       x.Output();
    588       DoLog(0) && (Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl);
    589       break;
    590 
    591     case 'c':
    592       DoLog(0) && (Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl);
    593       first = mol->AskAtom("Enter first atom: ");
    594       second = mol->AskAtom("Enter central atom: ");
    595       third  = mol->AskAtom("Enter last atom: ");
    596       tmp1 = tmp2 = tmp3 = 0.;
    597       x.CopyVector((const Vector *)&first->x);
    598       x.SubtractVector((const Vector *)&second->x);
    599       y.CopyVector((const Vector *)&third->x);
    600       y.SubtractVector((const Vector *)&second->x);
    601       DoLog(0) && (Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ");
    602       DoLog(0) && (Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl);
    603       break;
    604     case 'd':
    605       DoLog(0) && (Log() << Verbose(0) << "Evaluating prinicipal axis." << endl);
    606       DoLog(0) && (Log() << Verbose(0) << "Shall we rotate? [0/1]: ");
    607       cin >> Z;
    608       if ((Z >=0) && (Z <=1))
    609         mol->PrincipalAxisSystem((bool)Z);
    610       else
    611         mol->PrincipalAxisSystem(false);
    612       break;
    613     case 'e':
    614       {
    615         DoLog(0) && (Log() << Verbose(0) << "Evaluating volume of the convex envelope.");
    616         class Tesselation *TesselStruct = NULL;
    617         const LinkedCell *LCList = NULL;
    618         LCList = new LinkedCell(mol, 10.);
    619         FindConvexBorder(mol, TesselStruct, LCList, NULL);
    620         double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
    621         DoLog(0) && (Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl);\
    622         delete(LCList);
    623         delete(TesselStruct);
    624       }
    625       break;
    626     case 'f':
    627       mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps);
    628       break;
    629     case 'g':
    630       {
    631         char filename[255];
    632         DoLog(0) && (Log() << Verbose(0) << "Please enter filename: " << endl);
    633         cin >> filename;
    634         DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl);
    635         ofstream *output = new ofstream(filename, ios::trunc);
    636         if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
    637           DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl);
    638         else
    639           DoLog(2) && (Log() << Verbose(2) << "File stored." << endl);
    640         output->close();
    641         delete(output);
    642       }
    643       break;
    644   }
    645 };
    646 
    647 /** Submenu for measuring out the atoms in the molecule.
    648  * \param *mol molecule with all the atoms
    649  * \param *configuration configuration structure for the to be written config files of all fragments
    650  */
    651 static void FragmentAtoms(molecule *mol, config *configuration)
    652 {
    653   int Order1;
    654   clock_t start, end;
    655 
    656   DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl);
    657   DoLog(0) && (Log() << Verbose(0) << "What's the desired bond order: ");
    658   cin >> Order1;
    659   if (mol->first->next != mol->last) {  // there are bonds
    660     start = clock();
    661     mol->FragmentMolecule(Order1, configuration);
    662     end = clock();
    663     DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl);
    664   } else
    665     DoLog(0) && (Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl);
    666 };
    667 
    668 /********************************************** Submenu routine **************************************/
    669 
    670 /** Submenu for manipulating atoms.
    671  * \param *periode periodentafel
    672  * \param *molecules list of molecules whose atoms are to be manipulated
    673  */
    674 static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
    675 {
    676   atom *first, *second, *third;
    677   molecule *mol = NULL;
    678   Vector x,y,z,n; // coordinates for absolute point in cell volume
    679   double *factor; // unit factor if desired
    680   double bond, minBond;
    681   char choice;  // menu choice char
    682   bool valid;
    683 
    684   DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl);
    685   DoLog(0) && (Log() << Verbose(0) << "a - add an atom" << endl);
    686   DoLog(0) && (Log() << Verbose(0) << "r - remove an atom" << endl);
    687   DoLog(0) && (Log() << Verbose(0) << "b - scale a bond between atoms" << endl);
    688   DoLog(0) && (Log() << Verbose(0) << "t - turn an atom round another bond" << endl);
    689   DoLog(0) && (Log() << Verbose(0) << "u - change an atoms element" << endl);
    690   DoLog(0) && (Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl);
    691   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    692   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    693   if (molecules->NumberOfActiveMolecules() > 1)
    694     DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
    695   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    696   cin >> choice;
    697 
    698   switch (choice) {
    699     default:
    700       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    701       break;
    702 
    703     case 'a': // add atom
    704       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    705         if ((*ListRunner)->ActiveFlag) {
    706         mol = *ListRunner;
    707         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    708         AddAtoms(periode, mol);
    709       }
    710       break;
    711 
    712     case 'b': // scale a bond
    713       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    714         if ((*ListRunner)->ActiveFlag) {
    715         mol = *ListRunner;
    716         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    717         DoLog(0) && (Log() << Verbose(0) << "Scaling bond length between two atoms." << endl);
    718         first = mol->AskAtom("Enter first (fixed) atom: ");
    719         second = mol->AskAtom("Enter second (shifting) atom: ");
    720         minBond = 0.;
    721         for (int i=NDIM;i--;)
    722           minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
    723         minBond = sqrt(minBond);
    724         DoLog(0) && (Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl);
    725         DoLog(0) && (Log() << Verbose(0) << "Enter new bond length [a.u.]: ");
    726         cin >> bond;
    727         for (int i=NDIM;i--;) {
    728           second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
    729         }
    730         //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
    731         //second->Output(second->type->No, 1);
    732       }
    733       break;
    734 
    735     case 'c': // unit scaling of the metric
    736       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    737         if ((*ListRunner)->ActiveFlag) {
    738         mol = *ListRunner;
    739         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    740        DoLog(0) && (Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl);
    741        DoLog(0) && (Log() << Verbose(0) << "Enter three factors: ");
    742        factor = new double[NDIM];
    743        cin >> factor[0];
    744        cin >> factor[1];
    745        cin >> factor[2];
    746        valid = true;
    747        mol->Scale((const double ** const)&factor);
    748        delete[](factor);
    749       }
    750      break;
    751 
    752     case 'l': // measure distances or angles
    753       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    754         if ((*ListRunner)->ActiveFlag) {
    755         mol = *ListRunner;
    756         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    757         MeasureAtoms(periode, mol, configuration);
    758       }
    759       break;
    760 
    761     case 'r': // remove atom
    762       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    763         if ((*ListRunner)->ActiveFlag) {
    764         mol = *ListRunner;
    765         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    766         RemoveAtoms(mol);
    767       }
    768       break;
    769 
    770     case 't': // turn/rotate atom
    771       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    772         if ((*ListRunner)->ActiveFlag) {
    773           mol = *ListRunner;
    774           DoLog(0) && (Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl);
    775           first = mol->AskAtom("Enter turning atom: ");
    776           second = mol->AskAtom("Enter central atom: ");
    777           third  = mol->AskAtom("Enter bond atom: ");
    778           cout << Verbose(0) << "Enter new angle in degrees: ";
    779           double tmp = 0.;
    780           cin >> tmp;
    781           // calculate old angle
    782           x.CopyVector((const Vector *)&first->x);
    783           x.SubtractVector((const Vector *)&second->x);
    784           y.CopyVector((const Vector *)&third->x);
    785           y.SubtractVector((const Vector *)&second->x);
    786           double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
    787           cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
    788           cout << Verbose(0) << alpha << " degrees" << endl;
    789           // rotate
    790           z.MakeNormalVector(&x,&y);
    791           x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
    792           x.AddVector(&second->x);
    793           first->x.CopyVector(&x);
    794           // check new angle
    795           x.CopyVector((const Vector *)&first->x);
    796           x.SubtractVector((const Vector *)&second->x);
    797           alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
    798           cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
    799           cout << Verbose(0) << alpha << " degrees" << endl;
    800         }
    801       break;
    802 
    803     case 'u': // change an atom's element
    804       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    805         if ((*ListRunner)->ActiveFlag) {
    806         int Z;
    807         mol = *ListRunner;
    808         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    809         first = NULL;
    810         do {
    811           DoLog(0) && (Log() << Verbose(0) << "Change the element of which atom: ");
    812           cin >> Z;
    813         } while ((first = mol->FindAtom(Z)) == NULL);
    814         DoLog(0) && (Log() << Verbose(0) << "New element by atomic number Z: ");
    815         cin >> Z;
    816         first->type = periode->FindElement(Z);
    817         DoLog(0) && (Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl);
    818       }
    819       break;
    820   }
    821 };
    822 
    823 /** Submenu for manipulating molecules.
    824  * \param *periode periodentafel
    825  * \param *molecules list of molecule to manipulate
    826  */
    827 static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
    828 {
    829   atom *first = NULL;
    830   Vector x,y,z,n; // coordinates for absolute point in cell volume
    831   int j, axis, count, faktor;
    832   char choice;  // menu choice char
    833   molecule *mol = NULL;
    834   element **Elements;
    835   Vector **vectors;
    836   MoleculeLeafClass *Subgraphs = NULL;
    837 
    838   DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl);
    839   DoLog(0) && (Log() << Verbose(0) << "c - scale by unit transformation" << endl);
    840   DoLog(0) && (Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl);
    841   DoLog(0) && (Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl);
    842   DoLog(0) && (Log() << Verbose(0) << "g - center atoms in box" << endl);
    843   DoLog(0) && (Log() << Verbose(0) << "i - realign molecule" << endl);
    844   DoLog(0) && (Log() << Verbose(0) << "m - mirror all molecules" << endl);
    845   DoLog(0) && (Log() << Verbose(0) << "o - create connection matrix" << endl);
    846   DoLog(0) && (Log() << Verbose(0) << "t - translate molecule by vector" << endl);
    847   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    848   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    849   if (molecules->NumberOfActiveMolecules() > 1)
    850     DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
    851   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    852   cin >> choice;
    853 
    854   switch (choice) {
    855     default:
    856       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    857       break;
    858 
    859     case 'd': // duplicate the periodic cell along a given axis, given times
    860       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    861         if ((*ListRunner)->ActiveFlag) {
    862         mol = *ListRunner;
    863         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    864         DoLog(0) && (Log() << Verbose(0) << "State the axis [(+-)123]: ");
    865         cin >> axis;
    866         DoLog(0) && (Log() << Verbose(0) << "State the factor: ");
    867         cin >> faktor;
    868 
    869         mol->CountAtoms(); // recount atoms
    870         if (mol->getAtomCount() != 0) {  // if there is more than none
    871           count = mol->getAtomCount();  // is changed becausing of adding, thus has to be stored away beforehand
    872           Elements = new element *[count];
    873           vectors = new Vector *[count];
    874           j = 0;
    875           first = mol->start;
    876           while (first->next != mol->end) { // make a list of all atoms with coordinates and element
    877             first = first->next;
    878             Elements[j] = first->type;
    879             vectors[j] = &first->x;
    880             j++;
    881           }
    882           if (count != j)
    883             DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
    884           x.Zero();
    885           y.Zero();
    886           y.x[abs(axis)-1] = World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
    887           for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    888             x.AddVector(&y); // per factor one cell width further
    889             for (int k=count;k--;) { // go through every atom of the original cell
    890               first = new atom(); // create a new body
    891               first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    892               first->x.AddVector(&x);     // translate the coordinates
    893               first->type = Elements[k];  // insert original element
    894               mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
    895             }
    896           }
    897           if (mol->first->next != mol->last) // if connect matrix is present already, redo it
    898             mol->CreateAdjacencyList(mol->BondDistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
    899           // free memory
    900           delete[](Elements);
    901           delete[](vectors);
    902           // correct cell size
    903           if (axis < 0) { // if sign was negative, we have to translate everything
    904             x.Zero();
    905             x.AddVector(&y);
    906             x.Scale(-(faktor-1));
    907             mol->Translate(&x);
    908           }
    909           World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    910         }
    911       }
    912       break;
    913 
    914     case 'f':
    915       FragmentAtoms(mol, configuration);
    916       break;
    917 
    918     case 'g': // center the atoms
    919       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    920         if ((*ListRunner)->ActiveFlag) {
    921         mol = *ListRunner;
    922         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    923         CenterAtoms(mol);
    924       }
    925       break;
    926 
    927     case 'i': // align all atoms
    928       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    929         if ((*ListRunner)->ActiveFlag) {
    930         mol = *ListRunner;
    931         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    932         AlignAtoms(periode, mol);
    933       }
    934       break;
    935 
    936     case 'm': // mirror atoms along a given axis
    937       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    938         if ((*ListRunner)->ActiveFlag) {
    939         mol = *ListRunner;
    940         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    941         MirrorAtoms(mol);
    942       }
    943       break;
    944 
    945     case 'o': // create the connection matrix
    946       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    947         if ((*ListRunner)->ActiveFlag) {
    948           mol = *ListRunner;
    949           double bonddistance;
    950           clock_t start,end;
    951           DoLog(0) && (Log() << Verbose(0) << "What's the maximum bond distance: ");
    952           cin >> bonddistance;
    953           start = clock();
    954           mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
    955           end = clock();
    956           DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl);
    957         }
    958       break;
    959 
    960     case 't': // translate all atoms
    961       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    962         if ((*ListRunner)->ActiveFlag) {
    963         mol = *ListRunner;
    964         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    965         DoLog(0) && (Log() << Verbose(0) << "Enter translation vector." << endl);
    966         x.AskPosition(World::getInstance().getDomain(),0);
    967         mol->Center.AddVector((const Vector *)&x);
    968      }
    969      break;
    970   }
    971   // Free all
    972   if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
    973     while (Subgraphs->next != NULL) {
    974       Subgraphs = Subgraphs->next;
    975       delete(Subgraphs->previous);
    976     }
    977     delete(Subgraphs);
    978   }
    979 };
    980 
    981 
    982 /** Submenu for creating new molecules.
    983  * \param *periode periodentafel
    984  * \param *molecules list of molecules to add to
    985  */
    986 static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
    987 {
    988   char choice;  // menu choice char
    989   Vector center;
    990   int nr, count;
    991   molecule *mol = NULL;
    992 
    993   DoLog(0) && (Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl);
    994   DoLog(0) && (Log() << Verbose(0) << "c - create new molecule" << endl);
    995   DoLog(0) && (Log() << Verbose(0) << "l - load molecule from xyz file" << endl);
    996   DoLog(0) && (Log() << Verbose(0) << "n - change molecule's name" << endl);
    997   DoLog(0) && (Log() << Verbose(0) << "N - give molecules filename" << endl);
    998   DoLog(0) && (Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl);
    999   DoLog(0) && (Log() << Verbose(0) << "r - remove a molecule" << endl);
    1000   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    1001   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    1002   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    1003   cin >> choice;
    1004 
    1005   switch (choice) {
    1006     default:
    1007       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    1008       break;
    1009     case 'c':
    1010       mol = World::getInstance().createMolecule();
    1011       molecules->insert(mol);
    1012       break;
    1013 
    1014     case 'l': // load from XYZ file
    1015       {
    1016         char filename[MAXSTRINGSIZE];
    1017         DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl);
    1018         mol = World::getInstance().createMolecule();
    1019         do {
    1020           DoLog(0) && (Log() << Verbose(0) << "Enter file name: ");
    1021           cin >> filename;
    1022         } while (!mol->AddXYZFile(filename));
    1023         mol->SetNameFromFilename(filename);
    1024         // center at set box dimensions
    1025         mol->CenterEdge(&center);
    1026         double * const cell_size = World::getInstance().getDomain();
    1027         cell_size[0] = center.x[0];
    1028         cell_size[1] = 0;
    1029         cell_size[2] = center.x[1];
    1030         cell_size[3] = 0;
    1031         cell_size[4] = 0;
    1032         cell_size[5] = center.x[2];
    1033         molecules->insert(mol);
    1034       }
    1035       break;
    1036 
    1037     case 'n':
    1038       {
    1039         char filename[MAXSTRINGSIZE];
    1040         do {
    1041           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1042           cin >> nr;
    1043           mol = molecules->ReturnIndex(nr);
    1044         } while (mol == NULL);
    1045         DoLog(0) && (Log() << Verbose(0) << "Enter name: ");
    1046         cin >> filename;
    1047         strcpy(mol->name, filename);
    1048       }
    1049       break;
    1050 
    1051     case 'N':
    1052       {
    1053         char filename[MAXSTRINGSIZE];
    1054         do {
    1055           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1056           cin >> nr;
    1057           mol = molecules->ReturnIndex(nr);
    1058         } while (mol == NULL);
    1059         DoLog(0) && (Log() << Verbose(0) << "Enter name: ");
    1060         cin >> filename;
    1061         mol->SetNameFromFilename(filename);
    1062       }
    1063       break;
    1064 
    1065     case 'p': // parse XYZ file
    1066       {
    1067         char filename[MAXSTRINGSIZE];
    1068         mol = NULL;
    1069         do {
    1070           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1071           cin >> nr;
    1072           mol = molecules->ReturnIndex(nr);
    1073         } while (mol == NULL);
    1074         DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl);
    1075         do {
    1076           DoLog(0) && (Log() << Verbose(0) << "Enter file name: ");
    1077           cin >> filename;
    1078         } while (!mol->AddXYZFile(filename));
    1079         mol->SetNameFromFilename(filename);
    1080       }
    1081       break;
    1082 
    1083     case 'r':
    1084       DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1085       cin >> nr;
    1086       count = 1;
    1087       for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    1088         if (nr == (*ListRunner)->IndexNr) {
    1089           mol = *ListRunner;
    1090           molecules->ListOfMolecules.erase(ListRunner);
    1091           delete(mol);
    1092           break;
    1093         }
    1094       break;
    1095   }
    1096 };
    1097 
    1098 
    1099 /** Submenu for merging molecules.
    1100  * \param *periode periodentafel
    1101  * \param *molecules list of molecules to add to
    1102  */
    1103 static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
    1104 {
    1105   char choice;  // menu choice char
    1106 
    1107   DoLog(0) && (Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl);
    1108   DoLog(0) && (Log() << Verbose(0) << "a - simple add of one molecule to another" << endl);
    1109   DoLog(0) && (Log() << Verbose(0) << "b - count the number of bonds of two elements" << endl);
    1110   DoLog(0) && (Log() << Verbose(0) << "B - count the number of bonds of three elements " << endl);
    1111   DoLog(0) && (Log() << Verbose(0) << "e - embedding merge of two molecules" << endl);
    1112   DoLog(0) && (Log() << Verbose(0) << "h - count the number of hydrogen bonds" << endl);
    1113   DoLog(0) && (Log() << Verbose(0) << "b - count the number of hydrogen bonds" << endl);
    1114   DoLog(0) && (Log() << Verbose(0) << "m - multi-merge of all molecules" << endl);
    1115   DoLog(0) && (Log() << Verbose(0) << "s - scatter merge of two molecules" << endl);
    1116   DoLog(0) && (Log() << Verbose(0) << "t - simple merge of two molecules" << endl);
    1117   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    1118   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    1119   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    1120   cin >> choice;
    1121 
    1122   switch (choice) {
    1123     default:
    1124       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    1125       break;
    1126 
    1127     case 'a':
    1128       {
    1129         int src, dest;
    1130         molecule *srcmol = NULL, *destmol = NULL;
    1131         {
    1132           do {
    1133             DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: ");
    1134             cin >> dest;
    1135             destmol = molecules->ReturnIndex(dest);
    1136           } while ((destmol == NULL) && (dest != -1));
    1137           do {
    1138             DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to add from: ");
    1139             cin >> src;
    1140             srcmol = molecules->ReturnIndex(src);
    1141           } while ((srcmol == NULL) && (src != -1));
    1142           if ((src != -1) && (dest != -1))
    1143             molecules->SimpleAdd(srcmol, destmol);
    1144         }
    1145       }
    1146       break;
    1147 
    1148     case 'b':
    1149       {
    1150         const int nr = 2;
    1151         char *names[nr] = {"first", "second"};
    1152         int Z[nr];
    1153         element *elements[nr];
    1154         for (int i=0;i<nr;i++) {
    1155           Z[i] = 0;
    1156           do {
    1157             cout << "Enter " << names[i] << " element: ";
    1158             cin >> Z[i];
    1159           } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
    1160           elements[i] = periode->FindElement(Z[i]);
    1161         }
    1162         const int count = CountBondsOfTwo(molecules, elements[0], elements[1]);
    1163         cout << endl << "There are " << count << " ";
    1164         for (int i=0;i<nr;i++) {
    1165           if (i==0)
    1166             cout << elements[i]->symbol;
    1167           else
    1168             cout << "-" << elements[i]->symbol;
    1169         }
    1170         cout << " bonds." << endl;
    1171       }
    1172     break;
    1173 
    1174     case 'B':
    1175       {
    1176         const int nr = 3;
    1177         char *names[nr] = {"first", "second", "third"};
    1178         int Z[nr];
    1179         element *elements[nr];
    1180         for (int i=0;i<nr;i++) {
    1181           Z[i] = 0;
    1182           do {
    1183             cout << "Enter " << names[i] << " element: ";
    1184             cin >> Z[i];
    1185           } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
    1186           elements[i] = periode->FindElement(Z[i]);
    1187         }
    1188         const int count = CountBondsOfThree(molecules, elements[0], elements[1], elements[2]);
    1189         cout << endl << "There are " << count << " ";
    1190         for (int i=0;i<nr;i++) {
    1191           if (i==0)
    1192             cout << elements[i]->symbol;
    1193           else
    1194             cout << "-" << elements[i]->symbol;
    1195         }
    1196         cout << " bonds." << endl;
    1197       }
    1198     break;
    1199 
    1200     case 'e':
    1201       {
    1202         int src, dest;
    1203         molecule *srcmol = NULL, *destmol = NULL;
    1204         do {
    1205           DoLog(0) && (Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ");
    1206           cin >> src;
    1207           srcmol = molecules->ReturnIndex(src);
    1208         } while ((srcmol == NULL) && (src != -1));
    1209         do {
    1210           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ");
    1211           cin >> dest;
    1212           destmol = molecules->ReturnIndex(dest);
    1213         } while ((destmol == NULL) && (dest != -1));
    1214         if ((src != -1) && (dest != -1))
    1215           molecules->EmbedMerge(destmol, srcmol);
    1216       }
    1217       break;
    1218 
    1219     case 'h':
    1220       {
    1221         int Z;
    1222         cout << "Please enter interface element: ";
    1223         cin >> Z;
    1224         element * const InterfaceElement = periode->FindElement(Z);
    1225         cout << endl << "There are " << CountHydrogenBridgeBonds(molecules, InterfaceElement) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << "." << endl;
    1226       }
    1227       break;
    1228 
    1229     case 'm':
    1230       {
    1231         int nr;
    1232         molecule *mol = NULL;
    1233         do {
    1234           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into: ");
    1235           cin >> nr;
    1236           mol = molecules->ReturnIndex(nr);
    1237         } while ((mol == NULL) && (nr != -1));
    1238         if (nr != -1) {
    1239           int N = molecules->ListOfMolecules.size()-1;
    1240           int *src = new int(N);
    1241           for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    1242             if ((*ListRunner)->IndexNr != nr)
    1243               src[N++] = (*ListRunner)->IndexNr;       
    1244           molecules->SimpleMultiMerge(mol, src, N);
    1245           delete[](src);
    1246         }
    1247       }
    1248       break;
    1249 
    1250     case 's':
    1251       DoLog(0) && (Log() << Verbose(0) << "Not implemented yet." << endl);
    1252       break;
    1253 
    1254     case 't':
    1255       {
    1256         int src, dest;
    1257         molecule *srcmol = NULL, *destmol = NULL;
    1258         {
    1259           do {
    1260             DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: ");
    1261             cin >> dest;
    1262             destmol = molecules->ReturnIndex(dest);
    1263           } while ((destmol == NULL) && (dest != -1));
    1264           do {
    1265             DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to merge into: ");
    1266             cin >> src;
    1267             srcmol = molecules->ReturnIndex(src);
    1268           } while ((srcmol == NULL) && (src != -1));
    1269           if ((src != -1) && (dest != -1))
    1270             molecules->SimpleMerge(srcmol, destmol);
    1271         }
    1272       }
    1273       break;
    1274   }
    1275 };
    1276 
    1277 /********************************************** Test routine **************************************/
    1278 
    1279 /** Is called always as option 'T' in the menu.
    1280  * \param *molecules list of molecules
    1281  */
    1282 static void testroutine(MoleculeListClass *molecules)
    1283 {
    1284   // the current test routine checks the functionality of the KeySet&Graph concept:
    1285   // We want to have a multiindex (the KeySet) describing a unique subgraph
    1286   int i, comp, counter=0;
    1287 
    1288   // create a clone
    1289   molecule *mol = NULL;
    1290   if (molecules->ListOfMolecules.size() != 0) // clone
    1291     mol = (molecules->ListOfMolecules.front())->CopyMolecule();
    1292   else {
    1293     DoeLog(0) && (eLog()<< Verbose(0) << "I don't have anything to test on ... ");
    1294     performCriticalExit();
    1295     return;
    1296   }
    1297   atom *Walker = mol->start;
    1298 
    1299   // generate some KeySets
    1300   DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl);
    1301   KeySet TestSets[mol->getAtomCount()+1];
    1302   i=1;
    1303   while (Walker->next != mol->end) {
    1304     Walker = Walker->next;
    1305     for (int j=0;j<i;j++) {
    1306       TestSets[j].insert(Walker->nr);
    1307     }
    1308     i++;
    1309   }
    1310   DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl);
    1311   KeySetTestPair test;
    1312   test = TestSets[mol->getAtomCount()-1].insert(Walker->nr);
    1313   if (test.second) {
    1314     DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
    1315   } else {
    1316     DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl);
    1317   }
    1318   TestSets[mol->getAtomCount()].insert(mol->end->previous->nr);
    1319   TestSets[mol->getAtomCount()].insert(mol->end->previous->previous->previous->nr);
    1320 
    1321   // constructing Graph structure
    1322   DoLog(0) && (Log() << Verbose(0) << "Generating Subgraph class." << endl);
    1323   Graph Subgraphs;
    1324 
    1325   // insert KeySets into Subgraphs
    1326   DoLog(0) && (Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl);
    1327   for (int j=0;j<mol->getAtomCount();j++) {
    1328     Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
    1329   }
    1330   DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl);
    1331   GraphTestPair test2;
    1332   test2 = Subgraphs.insert(GraphPair (TestSets[mol->getAtomCount()],pair<int, double>(counter++, 1.)));
    1333   if (test2.second) {
    1334     DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
    1335   } else {
    1336     DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl);
    1337   }
    1338 
    1339   // show graphs
    1340   DoLog(0) && (Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl);
    1341   Graph::iterator A = Subgraphs.begin();
    1342   while (A !=  Subgraphs.end()) {
    1343     DoLog(0) && (Log() << Verbose(0) << (*A).second.first << ": ");
    1344     KeySet::iterator key = (*A).first.begin();
    1345     comp = -1;
    1346     while (key != (*A).first.end()) {
    1347       if ((*key) > comp)
    1348         DoLog(0) && (Log() << Verbose(0) << (*key) << " ");
    1349       else
    1350         DoLog(0) && (Log() << Verbose(0) << (*key) << "! ");
    1351       comp = (*key);
    1352       key++;
    1353     }
    1354     DoLog(0) && (Log() << Verbose(0) << endl);
    1355     A++;
    1356   }
    1357   delete(mol);
    1358 };
    1359 
    1360 
    1361 /** Tries given filename or standard on saving the config file.
    1362  * \param *ConfigFileName name of file
    1363  * \param *configuration pointer to configuration structure with all the values
    1364  * \param *periode pointer to periodentafel structure with all the elements
    1365  * \param *molecules list of molecules structure with all the atoms and coordinates
    1366  */
    1367 static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
    1368 {
    1369   char filename[MAXSTRINGSIZE];
    1370   ofstream output;
    1371   molecule *mol = World::getInstance().createMolecule();
    1372   mol->SetNameFromFilename(ConfigFileName);
    1373 
    1374   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    1375     DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
    1376   }
    1377 
    1378 
    1379   // first save as PDB data
    1380   if (ConfigFileName != NULL)
    1381     strcpy(filename, ConfigFileName);
    1382   if (output == NULL)
    1383     strcpy(filename,"main_pcp_linux");
    1384   DoLog(0) && (Log() << Verbose(0) << "Saving as pdb input ");
    1385   if (configuration->SavePDB(filename, molecules))
    1386     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    1387   else
    1388     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1389 
    1390   // then save as tremolo data file
    1391   if (ConfigFileName != NULL)
    1392     strcpy(filename, ConfigFileName);
    1393   if (output == NULL)
    1394     strcpy(filename,"main_pcp_linux");
    1395   DoLog(0) && (Log() << Verbose(0) << "Saving as tremolo data input ");
    1396   if (configuration->SaveTREMOLO(filename, molecules))
    1397     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    1398   else
    1399     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1400 
    1401   // translate each to its center and merge all molecules in MoleculeListClass into this molecule
    1402   int N = molecules->ListOfMolecules.size();
    1403   int *src = new int[N];
    1404   N=0;
    1405   for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
    1406     src[N++] = (*ListRunner)->IndexNr;
    1407     (*ListRunner)->Translate(&(*ListRunner)->Center);
    1408   }
    1409   molecules->SimpleMultiAdd(mol, src, N);
    1410   delete[](src);
    1411 
    1412   // ... and translate back
    1413   for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
    1414     (*ListRunner)->Center.Scale(-1.);
    1415     (*ListRunner)->Translate(&(*ListRunner)->Center);
    1416     (*ListRunner)->Center.Scale(-1.);
    1417   }
    1418 
    1419   DoLog(0) && (Log() << Verbose(0) << "Storing configuration ... " << endl);
    1420   // get correct valence orbitals
    1421   mol->CalculateOrbitals(*configuration);
    1422   configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
    1423   if (ConfigFileName != NULL) { // test the file name
    1424     strcpy(filename, ConfigFileName);
    1425     output.open(filename, ios::trunc);
    1426   } else if (strlen(configuration->configname) != 0) {
    1427     strcpy(filename, configuration->configname);
    1428     output.open(configuration->configname, ios::trunc);
    1429     } else {
    1430       strcpy(filename, DEFAULTCONFIG);
    1431       output.open(DEFAULTCONFIG, ios::trunc);
    1432     }
    1433   output.close();
    1434   output.clear();
    1435   DoLog(0) && (Log() << Verbose(0) << "Saving of config file ");
    1436   if (configuration->Save(filename, periode, mol))
    1437     DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
    1438   else
    1439     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1440 
    1441   // and save to xyz file
    1442   if (ConfigFileName != NULL) {
    1443     strcpy(filename, ConfigFileName);
    1444     strcat(filename, ".xyz");
    1445     output.open(filename, ios::trunc);
    1446   }
    1447   if (output == NULL) {
    1448     strcpy(filename,"main_pcp_linux");
    1449     strcat(filename, ".xyz");
    1450     output.open(filename, ios::trunc);
    1451   }
    1452   DoLog(0) && (Log() << Verbose(0) << "Saving of XYZ file ");
    1453   if (mol->MDSteps <= 1) {
    1454     if (mol->OutputXYZ(&output))
    1455       DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
    1456     else
    1457       DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1458   } else {
    1459     if (mol->OutputTrajectoriesXYZ(&output))
    1460       DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
    1461     else
    1462       DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1463   }
    1464   output.close();
    1465   output.clear();
    1466 
    1467   // and save as MPQC configuration
    1468   if (ConfigFileName != NULL)
    1469     strcpy(filename, ConfigFileName);
    1470   if (output == NULL)
    1471     strcpy(filename,"main_pcp_linux");
    1472   DoLog(0) && (Log() << Verbose(0) << "Saving as mpqc input ");
    1473   if (configuration->SaveMPQC(filename, mol))
    1474     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    1475   else
    1476     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1477 
    1478   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    1479     DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
    1480   }
    1481 
    1482   World::getInstance().destroyMolecule(mol);
    1483 };
    1484 
    1485 #endif
    1486 
    1487 /** Parses the command line options.
    1488  * Note that this function is from now on transitional. All commands that are not passed
    1489  * here are handled by CommandLineParser and the actions of CommandLineUIFactory.
    1490  * \param argc argument count
    1491  * \param **argv arguments array
    1492  * \param *molecules list of molecules structure
    1493  * \param *periode elements structure
    1494  * \param configuration config file structure
    1495  * \param *ConfigFileName pointer to config file name in **argv
    1496  * \param *PathToDatabases pointer to db's path in **argv
    1497  * \param &ArgcList list of arguments that we do not parse here
    1498  * \return exit code (0 - successful, all else - something's wrong)
    1499  */
    1500 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,
    1501                                    config& configuration, char **ConfigFileName, set<int> &ArgcList)
    1502 {
    1503   Vector x,y,z,n;  // coordinates for absolute point in cell volume
    1504   ifstream test;
    1505   ofstream output;
    1506   string line;
    1507   bool SaveFlag = false;
    1508   int ExitFlag = 0;
    1509   int j;
    1510   double volume = 0.;
    1511   enum ConfigStatus configPresent = absent;
    1512   int argptr;
    1513   molecule *mol = NULL;
    1514   string BondGraphFileName("\n");
    1515 
    1516   if (argc > 1) { // config file specified as option
    1517     // 1. : Parse options that just set variables or print help
    1518     argptr = 1;
    1519     do {
    1520       if (argv[argptr][0] == '-') {
    1521         DoLog(0) && (Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n");
    1522         argptr++;
    1523         switch(argv[argptr-1][1]) {
    1524           case 'h':
    1525           case 'H':
    1526           case '?':
    1527             ArgcList.insert(argptr-1);
    1528             return(1);
    1529             break;
    1530           case 'v':
    1531             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1532               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying verbosity: -v <level>" << endl);
    1533               performCriticalExit();
    1534             } else {
    1535               setVerbosity(atoi(argv[argptr]));
    1536               ArgcList.insert(argptr-1);
    1537               ArgcList.insert(argptr);
    1538               argptr++;
    1539             }
    1540             break;
    1541           case 'V':
    1542             ArgcList.insert(argptr-1);
    1543             return(1);
    1544             break;
    1545           case 'B':
    1546             if (ExitFlag == 0) ExitFlag = 1;
    1547             if ((argptr+5 >= argc)) {
    1548               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for setting Box: -B <xx> <<xy> <<xz> <yy> <yz> <zz>" << endl);
    1549               performCriticalExit();
    1550             } else {
    1551               ArgcList.insert(argptr-1);
    1552               ArgcList.insert(argptr);
    1553               ArgcList.insert(argptr+1);
    1554               ArgcList.insert(argptr+2);
    1555               ArgcList.insert(argptr+3);
    1556               ArgcList.insert(argptr+4);
    1557               ArgcList.insert(argptr+5);
    1558               argptr+=6;
    1559             }
    1560             break;
    1561           case 'e':
    1562             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1563               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl);
    1564               performCriticalExit();
    1565             } else {
    1566               ArgcList.insert(argptr-1);
    1567               ArgcList.insert(argptr);
    1568               argptr+=1;
    1569             }
    1570             break;
    1571           case 'g':
    1572             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1573               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl);
    1574               performCriticalExit();
    1575             } else {
    1576               ArgcList.insert(argptr-1);
    1577               ArgcList.insert(argptr);
    1578               argptr+=1;
    1579             }
    1580             break;
    1581           case 'M':
    1582             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1583               ExitFlag = 255;
    1584               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -M <basis name>" << endl);
    1585               performCriticalExit();
    1586             } else {
    1587               ArgcList.insert(argptr-1);
    1588               ArgcList.insert(argptr);
    1589               argptr+=1;
    1590             }
    1591             break;
    1592           case 'n':
    1593             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1594               ExitFlag = 255;
    1595               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting fast-parsing: -n <0/1>" << endl);
    1596               performCriticalExit();
    1597             } else {
    1598               ArgcList.insert(argptr-1);
    1599               ArgcList.insert(argptr);
    1600               argptr+=1;
    1601             }
    1602             break;
    1603           case 'X':
    1604             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1605               ExitFlag = 255;
    1606               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting default molecule name: -X <name>" << endl);
    1607               performCriticalExit();
    1608             } else {
    1609               ArgcList.insert(argptr-1);
    1610               ArgcList.insert(argptr);
    1611               argptr+=1;
    1612             }
    1613             break;
    1614           default:   // no match? Step on
    1615             argptr++;
    1616             break;
    1617         }
    1618       } else
    1619         argptr++;
    1620     } while (argptr < argc);
    1621 
    1622     // 3b. Find config file name and parse if possible, also BondGraphFileName
    1623     if (argv[1][0] != '-') {
    1624       // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
    1625       DoLog(0) && (Log() << Verbose(0) << "Config file given." << endl);
    1626       test.open(argv[1], ios::in);
    1627       if (test == NULL) {
    1628         //return (1);
    1629         output.open(argv[1], ios::out);
    1630         if (output == NULL) {
    1631           DoLog(1) && (Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl);
    1632           configPresent = absent;
    1633         } else {
    1634           DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl);
    1635           strcpy(*ConfigFileName, argv[1]);
    1636           configPresent = empty;
    1637           output.close();
    1638         }
    1639       } else {
    1640         test.close();
    1641         strcpy(*ConfigFileName, argv[1]);
    1642         DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... ");
    1643         switch (configuration.TestSyntax(*ConfigFileName, periode)) {
    1644           case 1:
    1645             DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl);
    1646             configuration.Load(*ConfigFileName, BondGraphFileName, periode, molecules);
    1647             configPresent = present;
    1648             break;
    1649           case 0:
    1650             DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl);
    1651             configuration.LoadOld(*ConfigFileName, BondGraphFileName, periode, molecules);
    1652             configPresent = present;
    1653             break;
    1654           default:
    1655             DoLog(0) && (Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl);
    1656             configPresent = empty;
    1657        }
    1658       }
    1659     } else
    1660       configPresent = absent;
    1661      // set mol to first active molecule
    1662      if (molecules->ListOfMolecules.size() != 0) {
    1663        for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    1664          if ((*ListRunner)->ActiveFlag) {
    1665            mol = *ListRunner;
    1666            break;
    1667          }
    1668      }
    1669      if (mol == NULL) {
    1670        mol = World::getInstance().createMolecule();
    1671        mol->ActiveFlag = true;
    1672        if (*ConfigFileName != NULL)
    1673          mol->SetNameFromFilename(*ConfigFileName);
    1674        molecules->insert(mol);
    1675      }
    1676 
    1677     // 4. parse again through options, now for those depending on elements db and config presence
    1678     argptr = 1;
    1679     do {
    1680       DoLog(0) && (Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl);
    1681       if (argv[argptr][0] == '-') {
    1682         argptr++;
    1683         if ((configPresent == present) || (configPresent == empty)) {
    1684           switch(argv[argptr-1][1]) {
    1685             case 'p':
    1686               if (ExitFlag == 0) ExitFlag = 1;
    1687               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1688                 ExitFlag = 255;
    1689                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl);
    1690                 performCriticalExit();
    1691               } else {
    1692                 SaveFlag = true;
    1693                 DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl);
    1694                 if (!mol->AddXYZFile(argv[argptr]))
    1695                   DoLog(2) && (Log() << Verbose(2) << "File not found." << endl);
    1696                 else {
    1697                   DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl);
    1698                   configPresent = present;
    1699                 }
    1700               }
    1701               break;
    1702             case 'a':
    1703               if (ExitFlag == 0) ExitFlag = 1;
    1704               if ((argptr+4 >= argc) || (argv[argptr][0] == '-')) {
    1705                 ExitFlag = 255;
    1706                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for adding atom: -a <Z> --position <x> <y> <z>" << endl);
    1707                 performCriticalExit();
    1708               } else {
    1709                 ArgcList.insert(argptr-1);
    1710                 ArgcList.insert(argptr);
    1711                 ArgcList.insert(argptr+1);
    1712                 ArgcList.insert(argptr+2);
    1713                 ArgcList.insert(argptr+3);
    1714                 ArgcList.insert(argptr+4);
    1715                 argptr+=5;
    1716               }
    1717               break;
    1718             default:   // no match? Don't step on (this is done in next switch's default)
    1719               break;
    1720           }
    1721         }
    1722         if (configPresent == present) {
    1723           switch(argv[argptr-1][1]) {
    1724             case 'D':
    1725               if (ExitFlag == 0) ExitFlag = 1;
    1726               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1727                 ExitFlag = 255;
    1728                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for depth-first-search analysis: -D <max. bond distance>" << endl);
    1729                 performCriticalExit();
    1730               } else {
    1731                 ArgcList.insert(argptr-1);
    1732                 ArgcList.insert(argptr);
    1733                 argptr+=1;
    1734               }
    1735               break;
    1736             case 'I':
    1737               DoLog(1) && (Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl);
    1738               ArgcList.insert(argptr-1);
    1739               argptr+=0;
    1740               break;
    1741             case 'C':
    1742               {
    1743                 if (ExitFlag == 0) ExitFlag = 1;
    1744                 if ((argptr+11 >= argc)) {
    1745                   ExitFlag = 255;
    1746                   DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C[p] <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl);
    1747                   performCriticalExit();
    1748                 } else {
    1749                   switch(argv[argptr][0]) {
    1750                     case 'E':
    1751                       ArgcList.insert(argptr-1);
    1752                       ArgcList.insert(argptr);
    1753                       ArgcList.insert(argptr+1);
    1754                       ArgcList.insert(argptr+2);
    1755                       ArgcList.insert(argptr+3);
    1756                       ArgcList.insert(argptr+4);
    1757                       ArgcList.insert(argptr+5);
    1758                       ArgcList.insert(argptr+6);
    1759                       ArgcList.insert(argptr+7);
    1760                       ArgcList.insert(argptr+8);
    1761                       ArgcList.insert(argptr+9);
    1762                       ArgcList.insert(argptr+10);
    1763                       ArgcList.insert(argptr+11);
    1764                       argptr+=12;
    1765                       break;
    1766 
    1767                     case 'P':
    1768                       ArgcList.insert(argptr-1);
    1769                       ArgcList.insert(argptr);
    1770                       ArgcList.insert(argptr+1);
    1771                       ArgcList.insert(argptr+2);
    1772                       ArgcList.insert(argptr+3);
    1773                       ArgcList.insert(argptr+4);
    1774                       ArgcList.insert(argptr+5);
    1775                       ArgcList.insert(argptr+6);
    1776                       ArgcList.insert(argptr+7);
    1777                       ArgcList.insert(argptr+8);
    1778                       ArgcList.insert(argptr+9);
    1779                       ArgcList.insert(argptr+10);
    1780                       ArgcList.insert(argptr+11);
    1781                       ArgcList.insert(argptr+12);
    1782                       ArgcList.insert(argptr+13);
    1783                       ArgcList.insert(argptr+14);
    1784                       argptr+=15;
    1785                       break;
    1786 
    1787                     case 'S':
    1788                       ArgcList.insert(argptr-1);
    1789                       ArgcList.insert(argptr);
    1790                       ArgcList.insert(argptr+1);
    1791                       ArgcList.insert(argptr+2);
    1792                       ArgcList.insert(argptr+3);
    1793                       ArgcList.insert(argptr+4);
    1794                       ArgcList.insert(argptr+5);
    1795                       ArgcList.insert(argptr+6);
    1796                       ArgcList.insert(argptr+7);
    1797                       ArgcList.insert(argptr+8);
    1798                       ArgcList.insert(argptr+9);
    1799                       ArgcList.insert(argptr+10);
    1800                       ArgcList.insert(argptr+11);
    1801                       ArgcList.insert(argptr+12);
    1802                       ArgcList.insert(argptr+13);
    1803                       ArgcList.insert(argptr+14);
    1804                       argptr+=15;
    1805                       break;
    1806 
    1807                     default:
    1808                       ExitFlag = 255;
    1809                       DoeLog(0) && (eLog()<< Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl);
    1810                       performCriticalExit();
    1811                       break;
    1812                   }
    1813                 }
    1814                 break;
    1815               }
    1816             case 'E':
    1817               if (ExitFlag == 0) ExitFlag = 1;
    1818               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr]))) {
    1819                 ExitFlag = 255;
    1820                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> --element <Z>" << endl);
    1821                 performCriticalExit();
    1822               } else {
    1823                 ArgcList.insert(argptr-1);
    1824                 ArgcList.insert(argptr);
    1825                 ArgcList.insert(argptr+1);
    1826                 ArgcList.insert(argptr+2);
    1827                 argptr+=3;
    1828               }
    1829               break;
    1830             case 'F':
    1831               if (ExitFlag == 0) ExitFlag = 1;
    1832               if ((argptr+12 >= argc) || (argv[argptr][0] == '-')) {
    1833                 ExitFlag = 255;
    1834                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling with molecule: -F <filler xyz file> --MaxDistance <distance or -1> --distances <x> <y> <z>  --lengths <surface> <randatm> <randmol> --DoRotate <0/1>" << endl);
    1835                 performCriticalExit();
    1836               } else {
    1837                 ArgcList.insert(argptr-1);
    1838                 ArgcList.insert(argptr);
    1839                 ArgcList.insert(argptr+1);
    1840                 ArgcList.insert(argptr+2);
    1841                 ArgcList.insert(argptr+3);
    1842                 ArgcList.insert(argptr+4);
    1843                 ArgcList.insert(argptr+5);
    1844                 ArgcList.insert(argptr+6);
    1845                 ArgcList.insert(argptr+7);
    1846                 ArgcList.insert(argptr+8);
    1847                 ArgcList.insert(argptr+9);
    1848                 ArgcList.insert(argptr+10);
    1849                 ArgcList.insert(argptr+11);
    1850                 ArgcList.insert(argptr+12);
    1851                 argptr+=13;
    1852               }
    1853               break;
    1854             case 'A':
    1855               if (ExitFlag == 0) ExitFlag = 1;
    1856               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1857                 ExitFlag =255;
    1858                 DoeLog(0) && (eLog()<< Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile> --molecule-by-id <molecule_id>" << endl);
    1859                 performCriticalExit();
    1860               } else {
    1861                 ArgcList.insert(argptr-1);
    1862                 ArgcList.insert(argptr);
    1863                 ArgcList.insert(argptr+1);
    1864                 ArgcList.insert(argptr+2);
    1865                 argptr+=3;
    1866               }
    1867               break;
    1868 
    1869             case 'J':
    1870               if (ExitFlag == 0) ExitFlag = 1;
    1871               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1872                 ExitFlag =255;
    1873                 DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of adjacency file: -J <path> --molecule-by-id <molecule_id>" << endl);
    1874                 performCriticalExit();
    1875               } else {
    1876                 ArgcList.insert(argptr-1);
    1877                 ArgcList.insert(argptr);
    1878                 ArgcList.insert(argptr+1);
    1879                 ArgcList.insert(argptr+2);
    1880                 argptr+=3;
    1881               }
    1882               break;
    1883 
    1884             case 'j':
    1885               if (ExitFlag == 0) ExitFlag = 1;
    1886               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1887                 ExitFlag =255;
    1888                 DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of bonds file: -j <path> --molecule-by-id <molecule_id>" << endl);
    1889                 performCriticalExit();
    1890               } else {
    1891                 ArgcList.insert(argptr-1);
    1892                 ArgcList.insert(argptr);
    1893                 ArgcList.insert(argptr+1);
    1894                 ArgcList.insert(argptr+2);
    1895                 argptr+=3;
    1896               }
    1897               break;
    1898 
    1899             case 'N':
    1900               if (ExitFlag == 0) ExitFlag = 1;
    1901               if ((argptr+4 >= argc) || (argv[argptr][0] == '-')){
    1902                 ExitFlag = 255;
    1903                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -N <molecule_id> --sphere-radius <radius> --nonconvex-file <output prefix>" << endl);
    1904                 performCriticalExit();
    1905               } else {
    1906                 ArgcList.insert(argptr-1);
    1907                 ArgcList.insert(argptr);
    1908                 ArgcList.insert(argptr+1);
    1909                 ArgcList.insert(argptr+2);
    1910                 ArgcList.insert(argptr+3);
    1911                 ArgcList.insert(argptr+4);
    1912                 argptr+=5;
    1913               }
    1914               break;
    1915             case 'S':
    1916               if (ExitFlag == 0) ExitFlag = 1;
    1917               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1918                 ExitFlag = 255;
    1919                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file> --molecule-by-id 0" << endl);
    1920                 performCriticalExit();
    1921               } else {
    1922                 ArgcList.insert(argptr-1);
    1923                 ArgcList.insert(argptr);
    1924                 ArgcList.insert(argptr+1);
    1925                 ArgcList.insert(argptr+2);
    1926                 argptr+=3;
    1927               }
    1928               break;
    1929             case 'L':
    1930               if (ExitFlag == 0) ExitFlag = 1;
    1931               if ((argptr+8 >= argc) || (argv[argptr][0] == '-')) {
    1932                 ExitFlag = 255;
    1933                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for linear interpolation: -L <prefix> --start-step <step0> --end-step <step1> --molecule-by-id 0 --id-mapping <0/1>" << endl);
    1934                 performCriticalExit();
    1935               } else {
    1936                 ArgcList.insert(argptr-1);
    1937                 ArgcList.insert(argptr);
    1938                 ArgcList.insert(argptr+1);
    1939                 ArgcList.insert(argptr+2);
    1940                 ArgcList.insert(argptr+3);
    1941                 ArgcList.insert(argptr+4);
    1942                 ArgcList.insert(argptr+5);
    1943                 ArgcList.insert(argptr+6);
    1944                 ArgcList.insert(argptr+7);
    1945                 ArgcList.insert(argptr+8);
    1946                 argptr+=9;
    1947               }
    1948               break;
    1949             case 'P':
    1950               if (ExitFlag == 0) ExitFlag = 1;
    1951               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1952                 ExitFlag = 255;
    1953                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file> --molecule-by-id <molecule_id>" << endl);
    1954                 performCriticalExit();
    1955               } else {
    1956                 ArgcList.insert(argptr-1);
    1957                 ArgcList.insert(argptr);
    1958                 ArgcList.insert(argptr+1);
    1959                 ArgcList.insert(argptr+2);
    1960                 argptr+=3;
    1961               }
    1962               break;
    1963             case 'R':
    1964               if (ExitFlag == 0) ExitFlag = 1;
    1965               if ((argptr+4 >= argc) || (argv[argptr][0] == '-'))  {
    1966                 ExitFlag = 255;
    1967                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <distance> --position <x> <y> <z>" << endl);
    1968                 performCriticalExit();
    1969               } else {
    1970                 ArgcList.insert(argptr-1);
    1971                 ArgcList.insert(argptr);
    1972                 ArgcList.insert(argptr+1);
    1973                 ArgcList.insert(argptr+2);
    1974                 ArgcList.insert(argptr+3);
    1975                 ArgcList.insert(argptr+4);
    1976                 argptr+=5;
    1977               }
    1978               break;
    1979             case 't':
    1980               if (ExitFlag == 0) ExitFlag = 1;
    1981               if ((argptr+4 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    1982                 ExitFlag = 255;
    1983                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z> --molecule-by-id <molecule_id> --periodic <0/1>" << endl);
    1984                 performCriticalExit();
    1985               } else {
    1986                 ArgcList.insert(argptr-1);
    1987                 ArgcList.insert(argptr);
    1988                 ArgcList.insert(argptr+1);
    1989                 ArgcList.insert(argptr+2);
    1990                 ArgcList.insert(argptr+3);
    1991                 ArgcList.insert(argptr+4);
    1992                 ArgcList.insert(argptr+5);
    1993                 ArgcList.insert(argptr+6);
    1994                 argptr+=7;
    1995               }
    1996               break;
    1997             case 's':
    1998               if (ExitFlag == 0) ExitFlag = 1;
    1999               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    2000                 ExitFlag = 255;
    2001                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl);
    2002                 performCriticalExit();
    2003               } else {
    2004                 ArgcList.insert(argptr-1);
    2005                 ArgcList.insert(argptr);
    2006                 ArgcList.insert(argptr+1);
    2007                 ArgcList.insert(argptr+2);
    2008                 argptr+=3;
    2009               }
    2010               break;
    2011             case 'b':
    2012               if (ExitFlag == 0) ExitFlag = 1;
    2013               if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
    2014                 ExitFlag = 255;
    2015                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
    2016                 performCriticalExit();
    2017               } else {
    2018                 ArgcList.insert(argptr-1);
    2019                 ArgcList.insert(argptr);
    2020                 ArgcList.insert(argptr+1);
    2021                 ArgcList.insert(argptr+2);
    2022                 ArgcList.insert(argptr+3);
    2023                 ArgcList.insert(argptr+4);
    2024                 ArgcList.insert(argptr+5);
    2025                 argptr+=6;
    2026               }
    2027               break;
    2028             case 'B':
    2029               if (ExitFlag == 0) ExitFlag = 1;
    2030               if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
    2031                 ExitFlag = 255;
    2032                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
    2033                 performCriticalExit();
    2034               } else {
    2035                 ArgcList.insert(argptr-1);
    2036                 ArgcList.insert(argptr);
    2037                 ArgcList.insert(argptr+1);
    2038                 ArgcList.insert(argptr+2);
    2039                 ArgcList.insert(argptr+3);
    2040                 ArgcList.insert(argptr+4);
    2041                 ArgcList.insert(argptr+5);
    2042                 argptr+=6;
    2043               }
    2044               break;
    2045             case 'c':
    2046               if (ExitFlag == 0) ExitFlag = 1;
    2047               if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    2048                 ExitFlag = 255;
    2049                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl);
    2050                 performCriticalExit();
    2051               } else {
    2052                 ArgcList.insert(argptr-1);
    2053                 ArgcList.insert(argptr);
    2054                 ArgcList.insert(argptr+1);
    2055                 ArgcList.insert(argptr+2);
    2056                 argptr+=3;
    2057               }
    2058               break;
    2059             case 'O':
    2060               if (ExitFlag == 0) ExitFlag = 1;
    2061               ArgcList.insert(argptr-1);
    2062               argptr+=0;
    2063               break;
    2064             case 'r':
    2065               if (ExitFlag == 0) ExitFlag = 1;
    2066               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])))  {
    2067                 ExitFlag = 255;
    2068                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl);
    2069                 performCriticalExit();
    2070               } else {
    2071                 ArgcList.insert(argptr-1);
    2072                 ArgcList.insert(argptr);
    2073                 argptr+=1;
    2074               }
    2075               break;
    2076             case 'f':
    2077               if (ExitFlag == 0) ExitFlag = 1;
    2078               if ((argptr+1 >= argc) || (argv[argptr][0] == '-')) {
    2079                 ExitFlag = 255;
    2080                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl);
    2081                 performCriticalExit();
    2082               } else {
    2083                 ArgcList.insert(argptr-1);
    2084                 ArgcList.insert(argptr);
    2085                 ArgcList.insert(argptr+1);
    2086                 ArgcList.insert(argptr+2);
    2087                 ArgcList.insert(argptr+3);
    2088                 ArgcList.insert(argptr+4);
    2089                 argptr+=5;
    2090               }
    2091               break;
    2092             case 'm':
    2093               if (ExitFlag == 0) ExitFlag = 1;
    2094               j = atoi(argv[argptr++]);
    2095               if ((j<0) || (j>1)) {
    2096                 DoeLog(1) && (eLog()<< Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl);
    2097                 j = 0;
    2098               }
    2099               if (j) {
    2100                 SaveFlag = true;
    2101                 DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl);
    2102                 mol->PrincipalAxisSystem((bool)j);
    2103               } else
    2104                 ArgcList.insert(argptr-1);
    2105                 argptr+=0;
    2106               break;
    2107             case 'o':
    2108               if (ExitFlag == 0) ExitFlag = 1;
    2109               if ((argptr+4 >= argc) || (argv[argptr][0] == '-')){
    2110                 ExitFlag = 255;
    2111                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <molecule_id> --output-file <output file> --output-file <binned output file>" << endl);
    2112                 performCriticalExit();
    2113               } else {
    2114                 ArgcList.insert(argptr-1);
    2115                 ArgcList.insert(argptr);
    2116                 ArgcList.insert(argptr+1);
    2117                 ArgcList.insert(argptr+2);
    2118                 ArgcList.insert(argptr+3);
    2119                 ArgcList.insert(argptr+4);
    2120                 argptr+=5;
    2121               }
    2122               break;
    2123             case 'U':
    2124               if (ExitFlag == 0) ExitFlag = 1;
    2125               if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
    2126                 ExitFlag = 255;
    2127                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl);
    2128                 performCriticalExit();
    2129               } else {
    2130                 volume = atof(argv[argptr++]);
    2131                 DoLog(0) && (Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl);
    2132               }
    2133             case 'u':
    2134               if (ExitFlag == 0) ExitFlag = 1;
    2135               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
    2136                 if (volume != -1)
    2137                   ExitFlag = 255;
    2138                   DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl);
    2139                   performCriticalExit();
    2140               } else {
    2141                 ArgcList.insert(argptr-1);
    2142                 ArgcList.insert(argptr);
    2143                 argptr+=1;
    2144               }
    2145               break;
    2146             case 'd':
    2147               if (ExitFlag == 0) ExitFlag = 1;
    2148               if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    2149                 ExitFlag = 255;
    2150                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl);
    2151                 performCriticalExit();
    2152               } else {
    2153                 ArgcList.insert(argptr-1);
    2154                 ArgcList.insert(argptr);
    2155                 ArgcList.insert(argptr+1);
    2156                 ArgcList.insert(argptr+2);
    2157                 argptr+=3;
    2158               }
    2159               break;
    2160             default:   // no match? Step on
    2161               if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
    2162                 argptr++;
    2163               break;
    2164           }
    2165         }
    2166       } else argptr++;
    2167     } while (argptr < argc);
    2168     if (SaveFlag)
    2169       configuration.SaveAll(*ConfigFileName, periode, molecules);
    2170   } else {  // no arguments, hence scan the elements db
    2171     if (periode->LoadPeriodentafel(configuration.databasepath))
    2172       DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl);
    2173     else
    2174       DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl);
    2175     configuration.RetrieveConfigPathAndName("main_pcp_linux");
    2176   }
    2177   return(ExitFlag);
    2178 };
    217980
    218081/********************************************** Main routine **************************************/
    218182
    218283void cleanUp(){
     84  FormatParserStorage::purgeInstance();
     85  ChangeTracker::purgeInstance();
    218386  World::purgeInstance();
    218487  logger::purgeInstance();
    … …  
    2197100int main(int argc, char **argv)
    2198101{
    2199     // while we are non interactive, we want to abort from asserts
    2200     //ASSERT_DO(Assert::Abort);
    2201     Vector x, y, z, n;
    2202     ifstream test;
    2203     ofstream output;
    2204     string line;
    2205     char **Arguments = NULL;
    2206     int ArgcSize = 0;
    2207     int ExitFlag = 0;
    2208     bool ArgumentsCopied = false;
    2209     char *ConfigFileName = new char[MAXSTRINGSIZE];
     102  // while we are non interactive, we want to abort from asserts
     103  //ASSERT_DO(Assert::Abort);
     104  string line;
     105  char **Arguments = NULL;
     106  int ArgcSize = 0;
     107  int ExitFlag = 0;
     108  bool ArgumentsCopied = false;
     109  std::string BondGraphFileName("\n");
     110  FormatParserStorage::getInstance().addMpqc();
     111  FormatParserStorage::getInstance().addPcp();
     112  FormatParserStorage::getInstance().addXyz();
    2210113
    2211     // print version check whether arguments are present at all
    2212     cout << ESPACKVersion << endl;
    2213     if (argc < 2) {
    2214       cout << "Obtain help with " << argv[0] << " -h." << endl;
    2215       cleanUp();
    2216       Memory::getState();
    2217       return(1);
    2218     }
     114  // print version check whether arguments are present at all
     115  cout << ESPACKVersion << endl;
     116  if (argc < 2) {
     117    cout << "Obtain help with " << argv[0] << " -h." << endl;
     118    cleanUp();
     119    Memory::getState();
     120    return(1);
     121  }
    2219122
    2220123
    2221     setVerbosity(0);
    2222     // need to init the history before any action is created
    2223     ActionHistory::init();
     124  setVerbosity(0);
     125  // need to init the history before any action is created
     126  ActionHistory::init();
    2224127
    2225     // In the interactive mode, we can leave the user the choice in case of error
    2226     ASSERT_DO(Assert::Ask);
     128  // In the interactive mode, we can leave the user the choice in case of error
     129  ASSERT_DO(Assert::Ask);
    2227130
    2228     // from this moment on, we need to be sure to deeinitialize in the correct order
    2229     // this is handled by the cleanup function
    2230     atexit(cleanUp);
     131  // from this moment on, we need to be sure to deeinitialize in the correct order
     132  // this is handled by the cleanup function
     133  atexit(cleanUp);
    2231134
    2232     // Parse command line options and if present create respective UI
    2233     {
    2234       set<int> ArgcList;
    2235       ArgcList.insert(0); // push back program!
    2236       ArgcList.insert(1); // push back config file name
    2237       // handle arguments by ParseCommandLineOptions()
    2238       ExitFlag = ParseCommandLineOptions(argc,argv,World::getInstance().getMolecules(),World::getInstance().getPeriode(),*World::getInstance().getConfig(), &ConfigFileName, ArgcList);
    2239       World::getInstance().setExitFlag(ExitFlag);
    2240       // copy all remaining arguments to a new argv
    2241       Arguments = new char *[ArgcList.size()];
    2242       cout << "The following arguments are handled by CommandLineParser: ";
    2243       for (set<int>::iterator ArgcRunner = ArgcList.begin(); ArgcRunner != ArgcList.end(); ++ArgcRunner) {
    2244         Arguments[ArgcSize] = new char[strlen(argv[*ArgcRunner])+2];
    2245         strcpy(Arguments[ArgcSize], argv[*ArgcRunner]);
    2246         cout << " " << argv[*ArgcRunner];
    2247         ArgcSize++;
    2248       }
    2249       cout << endl;
    2250       ArgumentsCopied = true;
    2251       // handle remaining arguments by CommandLineParser
    2252       MapOfActions::getInstance().AddOptionsToParser();
    2253       map <std::string, std::string> ShortFormToActionMap = MapOfActions::getInstance().getShortFormToActionMap();
    2254       CommandLineParser::getInstance().Run(ArgcSize,Arguments, ShortFormToActionMap);
    2255       if (!CommandLineParser::getInstance().isEmpty()) {
    2256         DoLog(0) && (Log() << Verbose(0) << "Setting UI to CommandLine." << endl);
    2257         UIFactory::registerFactory(new CommandLineUIFactory::description());
    2258         UIFactory::makeUserInterface("CommandLine");
     135  // Parse command line options and if present create respective UI
     136  {
     137    // construct bond graph
     138    if (World::getInstance().getConfig()->BG == NULL) {
     139      World::getInstance().getConfig()->BG = new BondGraph(World::getInstance().getConfig()->GetIsAngstroem());
     140      if (World::getInstance().getConfig()->BG->LoadBondLengthTable(BondGraphFileName)) {
     141        DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
    2259142      } else {
    2260         DoLog(0) && (Log() << Verbose(0) << "Setting UI to Text." << endl);
    2261         UIFactory::registerFactory(new TextUIFactory::description());
    2262         UIFactory::makeUserInterface("Text");
     143        DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
    2263144      }
    2264145    }
     146    // handle remaining arguments by CommandLineParser
     147    MapOfActions::getInstance().AddOptionsToParser();
     148    map <std::string, std::string> ShortFormToActionMap = MapOfActions::getInstance().getShortFormToActionMap();
     149    CommandLineParser::getInstance().Run(argc,argv, ShortFormToActionMap);
     150    if (!CommandLineParser::getInstance().isEmpty()) {
     151      DoLog(0) && (Log() << Verbose(0) << "Setting UI to CommandLine." << endl);
     152      UIFactory::registerFactory(new CommandLineUIFactory::description());
     153      UIFactory::makeUserInterface("CommandLine");
     154    } else {
     155      DoLog(0) && (Log() << Verbose(0) << "Setting UI to Text." << endl);
     156      UIFactory::registerFactory(new TextUIFactory::description());
     157      UIFactory::makeUserInterface("Text");
     158    }
     159  }
    2265160
    2266     {
    2267       MainWindow *mainWindow = UIFactory::getInstance().makeMainWindow();
    2268       mainWindow->display();
    2269       delete mainWindow;
    2270     }
     161  {
     162    MainWindow *mainWindow = UIFactory::getInstance().makeMainWindow();
     163    mainWindow->display();
     164    delete mainWindow;
     165  }
    2271166
    2272     Log() << Verbose(0) << "Saving to " << ConfigFileName << "." << endl;
    2273     World::getInstance().getConfig()->SaveAll(ConfigFileName, World::getInstance().getPeriode(), World::getInstance().getMolecules());
     167  FormatParserStorage::getInstance().SaveAll();
     168  ChangeTracker::getInstance().saveStatus();
    2274169
    2275170  // free the new argv
    … …  
    2279174    delete[](Arguments);
    2280175  }
    2281   delete[](ConfigFileName);
     176  //delete[](ConfigFileName);
    2282177
    2283178  ExitFlag = World::getInstance().getExitFlag();
Note: See TracChangeset for help on using the changeset viewer.