Changeset a356f2


Ignore:
Timestamp:
Jun 25, 2010, 8:31:18 AM (15 years ago)
Author:
Frederik Heber <heber@…>
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:
0d1ad0, 326bbe
Parents:
4f7f34e
Message:

Cleaned up builder.cpp

  • ParseComandLineOptions() along with the already #ifdef 0'd functions removed from builder.cpp
  • DOCU: Rewrote introductory to Molecuilder.

Signed-off-by: Frederik Heber <heber@…>

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

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