Changeset 04488a for src/builder.cpp


Ignore:
Timestamp:
Jun 25, 2010, 9:57:15 AM (15 years ago)
Author:
Tillmann Crueger <crueger@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
b6dbff
Parents:
ce4487 (diff), 0d1ad0 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'StructureRefactoring' into QT4Refactoring

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/builder.cpp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

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