Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

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