source: src/builder.cpp@ 65b6e0

/** \file builder.cpp
 *
 * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
 * The output is the complete configuration file for PCP for direct use.
 * Features:
 * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
 * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
 *
 */

/*! \mainpage Molecuilder - a molecular set builder
 *
 * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
 *
 * \section about About the Program
 *
 *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
 *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
 *  already constructed atoms.
 *
 *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
 *  molecular dynamics implementation.
 *
 * \section install Installation
 *
 *  Installation should without problems succeed as follows:
 *  -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
 *  -# make
 *  -# make install
 *
 *  Further useful commands are
 *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
 *  -# make doxygen-doc: Creates these html pages out of the documented source
 *
 * \section run Running
 *
 *  The program can be executed by running: ./molecuilder
 *
 *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
 *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
 *  later re-execution.
 *
 * \section ref References
 *
 *  For the special configuration file format, see the documentation of pcp.
 *
 */


using namespace std;

#include "analysis_correlation.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "bondgraph.hpp"
#include "boundary.hpp"
#include "config.hpp"
#include "element.hpp"
#include "ellipsoid.hpp"
#include "helpers.hpp"
#include "leastsquaremin.hpp"
#include "linkedcell.hpp"
#include "log.hpp"
#include "memoryusageobserverunittest.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "menu.hpp"

/** Parses the command line options.
 * \param argc argument count
 * \param **argv arguments array
 * \param *molecules list of molecules structure
 * \param *periode elements structure
 * \param configuration config file structure
 * \param *ConfigFileName pointer to config file name in **argv
 * \param *PathToDatabases pointer to db's path in **argv
 * \return exit code (0 - successful, all else - something's wrong)
 */
static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,\
                                   config& configuration, char *&ConfigFileName, oldmenu *main_menu)
{
  Vector x,y,z,n;  // coordinates for absolute point in cell volume
  double *factor; // unit factor if desired
  ifstream test;
  ofstream output;
  string line;
  atom *first;
  bool SaveFlag = false;
  int ExitFlag = 0;
  int j;
  double volume = 0.;
  enum ConfigStatus configPresent = absent;
  clock_t start,end;
  int argptr;
  molecule *mol = NULL;
  string BondGraphFileName("");
  int verbosity = 0;
  strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);

  if (argc > 1) { // config file specified as option
    // 1. : Parse options that just set variables or print help
    argptr = 1;
    do {
      if (argv[argptr][0] == '-') {
        Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
        argptr++;
        switch(argv[argptr-1][1]) {
          case 'h':
          case 'H':
          case '?':
            Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl;
            Log() << Verbose(0) << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
            Log() << Verbose(0) << "or simply " << argv[0] << " without arguments for interactive session." << endl;
            Log() << Verbose(0) << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
            Log() << Verbose(0) << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
            Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
            Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
            Log() << Verbose(0) << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
            Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
            Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
            Log() << Verbose(0) << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl;
            Log() << Verbose(0) << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
            Log() << Verbose(0) << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
            Log() << Verbose(0) << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
            Log() << Verbose(0) << "\t-O\tCenter atoms in origin." << endl;
            Log() << Verbose(0) << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
            Log() << Verbose(0) << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
            Log() << Verbose(0) << "\t-r <id>\t\tRemove an atom with given id." << endl;
            Log() << Verbose(0) << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl;
            Log() << Verbose(0) << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-S <file> Store temperatures from the config file in <file>." << endl;
            Log() << Verbose(0) << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl;
            Log() << Verbose(0) << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
            Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl;
            Log() << Verbose(0) << "\t-V\t\tGives version information." << endl;
            Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl;
            return (1);
            break;
          case 'v':
            while (argv[argptr-1][verbosity+1] == 'v') {
              verbosity++;
            }
            setVerbosity(verbosity);
            Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl;
            break;
          case 'V':
            Log() << Verbose(0) << argv[0] << " " << VERSIONSTRING << endl;
            Log() << Verbose(0) << "Build your own molecule position set." << endl;
            return (1);
            break;
          case 'e':
            if ((argptr >= argc) || (argv[argptr][0] == '-')) {
              eLog() << Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
              performCriticalExit();
            } else {
              Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl;
              strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1);
              argptr+=1;
            }
            break;
          case 'g':
            if ((argptr >= argc) || (argv[argptr][0] == '-')) {
              eLog() << Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl;
              performCriticalExit();
            } else {
              BondGraphFileName = argv[argptr];
              Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl;
              argptr+=1;
            }
            break;
          case 'n':
            Log() << Verbose(0) << "I won't parse trajectories." << endl;
            configuration.FastParsing = true;
            break;
          default:   // no match? Step on
            argptr++;
            break;
        }
      } else
        argptr++;
    } while (argptr < argc);

    // 3a. Parse the element database
    if (periode->LoadPeriodentafel(configuration.databasepath)) {
      Log() << Verbose(0) << "Element list loaded successfully." << endl;
      //periode->Output();
    } else {
      Log() << Verbose(0) << "Element list loading failed." << endl;
      return 1;
    }
    // 3b. Find config file name and parse if possible, also BondGraphFileName
    if (argv[1][0] != '-') {
      // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
      Log() << Verbose(0) << "Config file given." << endl;
      test.open(argv[1], ios::in);
      if (test == NULL) {
        //return (1);
        output.open(argv[1], ios::out);
        if (output == NULL) {
          Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
          configPresent = absent;
        } else {
          Log() << Verbose(0) << "Empty configuration file." << endl;
          ConfigFileName = argv[1];
          configPresent = empty;
          output.close();
        }
      } else {
        test.close();
        ConfigFileName = argv[1];
        Log() << Verbose(1) << "Specified config file found, parsing ... ";
        switch (configuration.TestSyntax(ConfigFileName, periode)) {
          case 1:
            Log() << Verbose(0) << "new syntax." << endl;
            configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules);
            configPresent = present;
            break;
          case 0:
            Log() << Verbose(0) << "old syntax." << endl;
            configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules);
            configPresent = present;
            break;
          default:
            Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl;
            configPresent = empty;
       }
      }
    } else
      configPresent = absent;
     // set mol to first active molecule
     if (molecules->ListOfMolecules.size() != 0) {
       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
         if ((*ListRunner)->ActiveFlag) {
           mol = *ListRunner;
           break;
         }
     }
     if (mol == NULL) {
       mol = new molecule(periode);
       mol->ActiveFlag = true;
       molecules->insert(mol);
     }

    // 4. parse again through options, now for those depending on elements db and config presence
    argptr = 1;
    do {
      Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl;
      if (argv[argptr][0] == '-') {
        argptr++;
        if ((configPresent == present) || (configPresent == empty)) {
          switch(argv[argptr-1][1]) {
            case 'p':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl;
                if (!mol->AddXYZFile(argv[argptr]))
                  Log() << Verbose(2) << "File not found." << endl;
                else {
                  Log() << Verbose(2) << "File found and parsed." << endl;
                  configPresent = present;
                }
              }
              break;
            case 'a':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3]))) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
                first = new atom;
                first->type = periode->FindElement(atoi(argv[argptr]));
                if (first->type != NULL)
                  Log() << Verbose(2) << "found element " << first->type->name << endl;
                for (int i=NDIM;i--;)
                  first->x.x[i] = atof(argv[argptr+1+i]);
                if (first->type != NULL) {
                  mol->AddAtom(first);  // add to molecule
                  if ((configPresent == empty) && (mol->AtomCount != 0))
                    configPresent = present;
                } else
                  eLog() << Verbose(1) << "Could not find the specified element." << endl;
                argptr+=4;
              }
              break;
            default:   // no match? Don't step on (this is done in next switch's default)
              break;
          }
        }
        if (configPresent == present) {
          switch(argv[argptr-1][1]) {
            case 'M':
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;
                performCriticalExit();
              } else {
                configuration.basis = argv[argptr];
                Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;
                argptr+=1;
              }
              break;
            case 'D':
              if (ExitFlag == 0) ExitFlag = 1;
              {
                Log() << Verbose(1) << "Depth-First-Search Analysis." << endl;
                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
                int *MinimumRingSize = new int[mol->AtomCount];
                atom ***ListOfLocalAtoms = NULL;
                class StackClass<bond *> *BackEdgeStack = NULL;
                class StackClass<bond *> *LocalBackEdgeStack = NULL;
                mol->CreateAdjacencyList(atof(argv[argptr]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
                Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
                if (Subgraphs != NULL) {
                  int FragmentCounter = 0;
                  while (Subgraphs->next != NULL) {
                    Subgraphs = Subgraphs->next;
                    Subgraphs->FillBondStructureFromReference(mol, FragmentCounter, ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
                    LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
                    Subgraphs->Leaf->PickLocalBackEdges(ListOfLocalAtoms[FragmentCounter], BackEdgeStack, LocalBackEdgeStack);
                    Subgraphs->Leaf->CyclicStructureAnalysis(LocalBackEdgeStack, MinimumRingSize);
                    delete(LocalBackEdgeStack);
                    delete(Subgraphs->previous);
                    FragmentCounter++;
                  }
                  delete(Subgraphs);
                  for (int i=0;i<FragmentCounter;i++)
                    Free(&ListOfLocalAtoms[i]);
                  Free(&ListOfLocalAtoms);
                }
                delete(BackEdgeStack);
                delete[](MinimumRingSize);
              }
              //argptr+=1;
              break;
            case 'C':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = false;
                ofstream output(argv[argptr+1]);
                ofstream binoutput(argv[argptr+2]);
                const double radius = 5.;

                // get the boundary
                class molecule *Boundary = NULL;
                class Tesselation *TesselStruct = NULL;
                const LinkedCell *LCList = NULL;
                // find biggest molecule
                int counter  = 0;
                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                  if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
                    Boundary = *BigFinder;
                  }
                  counter++;
                }
                bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
                counter = 0;
                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                  Actives[counter] = (*BigFinder)->ActiveFlag;
                  (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
                }
                LCList = new LinkedCell(Boundary, 2.*radius);
                element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
                FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
                int ranges[NDIM] = {1,1,1};
                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
                OutputCorrelation ( &binoutput, binmap );
                output.close();
                binoutput.close();
                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
                  (*BigFinder)->ActiveFlag = Actives[counter];
                Free(&Actives);
                delete(LCList);
                delete(TesselStruct);
                argptr+=3;
              }
              break;
            case 'E':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
                first = mol->FindAtom(atoi(argv[argptr]));
                first->type = periode->FindElement(atoi(argv[argptr+1]));
                argptr+=2;
              }
              break;
            case 'F':
              if (ExitFlag == 0) ExitFlag = 1;
              if (argptr+5 >=argc) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Filling Box with water molecules." << endl;
                // construct water molecule
                molecule *filler = new molecule(periode);;
                molecule *Filling = NULL;
                atom *second = NULL, *third = NULL;
                first = new atom();
                first->type = periode->FindElement(1);
                first->x.Init(0.441, -0.143, 0.);
                filler->AddAtom(first);
                second = new atom();
                second->type = periode->FindElement(1);
                second->x.Init(-0.464, 1.137, 0.0);
                filler->AddAtom(second);
                third = new atom();
                third->type = periode->FindElement(8);
                third->x.Init(-0.464, 0.177, 0.);
                filler->AddAtom(third);
                filler->AddBond(first, third, 1);
                filler->AddBond(second, third, 1);
                // call routine
                double distance[NDIM];
                for (int i=0;i<NDIM;i++)
                  distance[i] = atof(argv[argptr+i]);
                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
                if (Filling != NULL) {
                  molecules->insert(Filling);
                }
                delete(filler);
                argptr+=6;
              }
              break;
            case 'A':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag =255;
                eLog() << Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
                performCriticalExit();
              } else {
                Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl;
                ifstream *input = new ifstream(argv[argptr]);
                mol->CreateAdjacencyListFromDbondFile(input);
                input->close();
                argptr+=1;
              }
              break;
            case 'N':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
                performCriticalExit();
              } else {
                class Tesselation *T = NULL;
                const LinkedCell *LCList = NULL;
                molecule * Boundary = NULL;
                //string filename(argv[argptr+1]);
                //filename.append(".csv");
                Log() << Verbose(0) << "Evaluating non-convex envelope of biggest molecule.";
                Log() << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
                // find biggest molecule
                int counter  = 0;
                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
                  (*BigFinder)->CountAtoms();
                  if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
                    Boundary = *BigFinder;
                  }
                  counter++;
                }
                Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl;
                start = clock();
                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
                FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]);
                //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
                end = clock();
                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
                delete(LCList);
                argptr+=2;
              }
              break;
            case 'S':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl;
                performCriticalExit();
              } else {
                Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
                ofstream *output = new ofstream(argv[argptr], ios::trunc);
                if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
                  Log() << Verbose(2) << "File could not be written." << endl;
                else
                  Log() << Verbose(2) << "File stored." << endl;
                output->close();
                delete(output);
                argptr+=1;
              }
              break;
            case 'L':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
                if (atoi(argv[argptr+3]) == 1)
                  Log() << Verbose(1) << "Using Identity for the permutation map." << endl;
                if (!mol->LinearInterpolationBetweenConfiguration(atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false)
                  Log() << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
                else
                  Log() << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
                argptr+=4;
              }
              break;
            case 'P':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
                if (!mol->VerletForceIntegration(argv[argptr], configuration))
                  Log() << Verbose(2) << "File not found." << endl;
                else
                  Log() << Verbose(2) << "File found and parsed." << endl;
                argptr+=1;
              }
              break;
            case 'R':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl;
                double tmp1 = atof(argv[argptr+1]);
                atom *third = mol->FindAtom(atoi(argv[argptr]));
                atom *first = mol->start;
                if ((third != NULL) && (first != mol->end)) {
                  atom *second = first->next;
                  while(second != mol->end) {
                    first = second;
                    second = first->next;
                    if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
                      mol->RemoveAtom(first);
                  }
                } else {
                  eLog() << Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl;
                }
                argptr+=2;
              }
              break;
            case 't':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
                performCriticalExit();
              } else {
                if (ExitFlag == 0) ExitFlag = 1;
                SaveFlag = true;
                Log() << Verbose(1) << "Translating all ions by given vector." << endl;
                for (int i=NDIM;i--;)
                  x.x[i] = atof(argv[argptr+i]);
                mol->Translate((const Vector *)&x);
                argptr+=3;
              }
              break;
            case 'T':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl;
                performCriticalExit();
              } else {
                if (ExitFlag == 0) ExitFlag = 1;
                SaveFlag = true;
                Log() << Verbose(1) << "Translating all ions periodically by given vector." << endl;
                for (int i=NDIM;i--;)
                  x.x[i] = atof(argv[argptr+i]);
                mol->TranslatePeriodically((const Vector *)&x);
                argptr+=3;
              }
              break;
            case 's':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Scaling all ion positions by factor." << endl;
                factor = new double[NDIM];
                factor[0] = atof(argv[argptr]);
                factor[1] = atof(argv[argptr+1]);
                factor[2] = atof(argv[argptr+2]);
                mol->Scale((const double ** const)&factor);
                for (int i=0;i<NDIM;i++) {
                  j += i+1;
                  x.x[i] = atof(argv[NDIM+i]);
                  mol->cell_size[j]*=factor[i];
                }
                delete[](factor);
                argptr+=3;
              }
              break;
            case 'b':
              if (ExitFlag == 0) ExitFlag = 1;
              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])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
                for (int i=0;i<6;i++) {
                  mol->cell_size[i] = atof(argv[argptr+i]);
                }
                // center
                mol->CenterInBox();
                argptr+=6;
              }
              break;
            case 'B':
              if (ExitFlag == 0) ExitFlag = 1;
              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])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
                for (int i=0;i<6;i++) {
                  mol->cell_size[i] = atof(argv[argptr+i]);
                }
                // center
                mol->BoundInBox();
                argptr+=6;
              }
              break;
            case 'c':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                j = -1;
                Log() << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
                // make every coordinate positive
                mol->CenterEdge(&x);
                // update Box of atoms by boundary
                mol->SetBoxDimension(&x);
                // translate each coordinate by boundary
                j=-1;
                for (int i=0;i<NDIM;i++) {
                  j += i+1;
                  x.x[i] = atof(argv[argptr+i]);
                  mol->cell_size[j] += x.x[i]*2.;
                }
                mol->Translate((const Vector *)&x);
                argptr+=3;
              }
              break;
            case 'O':
              if (ExitFlag == 0) ExitFlag = 1;
              SaveFlag = true;
              Log() << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl;
              x.Zero();
              mol->CenterEdge(&x);
              mol->SetBoxDimension(&x);
              argptr+=0;
              break;
            case 'r':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])))  {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl;
                atom *first = mol->FindAtom(atoi(argv[argptr]));
                mol->RemoveAtom(first);
                argptr+=1;
              }
              break;
            case 'f':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
                performCriticalExit();
              } else {
                Log() << Verbose(0) << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
                Log() << Verbose(0) << "Creating connection matrix..." << endl;
                start = clock();
                mol->CreateAdjacencyList(atof(argv[argptr++]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
                Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
                if (mol->first->next != mol->last) {
                  ExitFlag = mol->FragmentMolecule(atoi(argv[argptr]), &configuration);
                }
                end = clock();
                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
                argptr+=2;
              }
              break;
            case 'm':
              if (ExitFlag == 0) ExitFlag = 1;
              j = atoi(argv[argptr++]);
              if ((j<0) || (j>1)) {
                eLog() << Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
                j = 0;
              }
              if (j) {
                SaveFlag = true;
                Log() << Verbose(0) << "Converting to prinicipal axis system." << endl;
              } else
                Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;
              mol->PrincipalAxisSystem((bool)j);
              break;
            case 'o':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;
                performCriticalExit();
              } else {
                class Tesselation *TesselStruct = NULL;
                const LinkedCell *LCList = NULL;
                Log() << Verbose(0) << "Evaluating volume of the convex envelope.";
                Log() << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;
                Log() << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;
                LCList = new LinkedCell(mol, 10.);
                //FindConvexBorder(mol, LCList, argv[argptr]);
                FindNonConvexBorder(mol, TesselStruct, LCList, 5., argv[argptr+1]);
//                RemoveAllBoundaryPoints(TesselStruct, mol, argv[argptr]);
                double volumedifference = ConvexizeNonconvexEnvelope(TesselStruct, mol, argv[argptr]);
                double clustervolume = VolumeOfConvexEnvelope(TesselStruct, &configuration);
                Log() << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                Log() << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                delete(TesselStruct);
                delete(LCList);
                argptr+=2;
              }
              break;
            case 'U':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
                performCriticalExit();
              } else {
                volume = atof(argv[argptr++]);
                Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
              }
            case 'u':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
                if (volume != -1)
                  ExitFlag = 255;
                  eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
                  performCriticalExit();
              } else {
                double density;
                SaveFlag = true;
                Log() << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
                density = atof(argv[argptr++]);
                if (density < 1.0) {
                  eLog() << Verbose(1) << "Density must be greater than 1.0g/cm^3 !" << endl;
                  density = 1.3;
                }
//                for(int i=0;i<NDIM;i++) {
//                  repetition[i] = atoi(argv[argptr++]);
//                  if (repetition[i] < 1)
//                    eLog() << Verbose(1) << "repetition value must be greater 1!" << endl;
//                  repetition[i] = 1;
//                }
                PrepareClustersinWater(&configuration, mol, volume, density);  // if volume == 0, will calculate from ConvexEnvelope
              }
              break;
            case 'd':
              if (ExitFlag == 0) ExitFlag = 1;
              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                ExitFlag = 255;
                eLog() << Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
                performCriticalExit();
              } else {
                SaveFlag = true;
                for (int axis = 1; axis <= NDIM; axis++) {
                  int faktor = atoi(argv[argptr++]);
                  int count;
                  element ** Elements;
                  Vector ** vectors;
                  if (faktor < 1) {
                    eLog() << Verbose(1) << "Repetition factor mus be greater than 1!" << endl;
                    faktor = 1;
                  }
                  mol->CountAtoms();  // recount atoms
                  if (mol->AtomCount != 0) {  // if there is more than none
                    count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
                    Elements = new element *[count];
                    vectors = new Vector *[count];
                    j = 0;
                    first = mol->start;
                    while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
                      first = first->next;
                      Elements[j] = first->type;
                      vectors[j] = &first->x;
                      j++;
                    }
                    if (count != j)
                      eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
                    x.Zero();
                    y.Zero();
                    y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
                    for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
                      x.AddVector(&y); // per factor one cell width further
                      for (int k=count;k--;) { // go through every atom of the original cell
                        first = new atom(); // create a new body
                        first->x.CopyVector(vectors[k]);  // use coordinate of original atom
                        first->x.AddVector(&x);      // translate the coordinates
                        first->type = Elements[k];  // insert original element
                        mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
                      }
                    }
                    // free memory
                    delete[](Elements);
                    delete[](vectors);
                    // correct cell size
                    if (axis < 0) { // if sign was negative, we have to translate everything
                      x.Zero();
                      x.AddVector(&y);
                      x.Scale(-(faktor-1));
                      mol->Translate(&x);
                    }
                    mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
                  }
                }
              }
              break;
            default:   // no match? Step on
              if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
                argptr++;
              break;
          }
        }
      } else argptr++;
    } while (argptr < argc);
    if (SaveFlag)
      main_menu->SaveConfig(ConfigFileName, &configuration, periode, molecules);
  } else {  // no arguments, hence scan the elements db
    if (periode->LoadPeriodentafel(configuration.databasepath))
      Log() << Verbose(0) << "Element list loaded successfully." << endl;
    else
      Log() << Verbose(0) << "Element list loading failed." << endl;
    configuration.RetrieveConfigPathAndName("main_pcp_linux");
  }
  return(ExitFlag);
};

/********************************************** Main routine **************************************/

int main(int argc, char **argv)
{
  periodentafel *periode = new periodentafel;
    MoleculeListClass *molecules = new MoleculeListClass;
    molecule *mol = NULL;
    config *configuration = new config;
    Vector x, y, z, n;
    ifstream test;
    ofstream output;
    string line;
    char *ConfigFileName = NULL;
    int j;
    oldmenu *main_menu;
    main_menu = new oldmenu;
    setVerbosity(0);
    /* main menu is needed to call actions inside the menu */
    /* structure of ParseCommandLineOptions will be refactored later */
    j = ParseCommandLineOptions(argc, argv, molecules, periode, *configuration, ConfigFileName,main_menu);
    switch (j){
        case 255:
        case 2:
        case 1:
            delete (molecules);
            delete (periode);
            delete (configuration);
            Log() << Verbose(0) << "Maximum of allocated memory: " << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
            Log() << Verbose(0) << "Remaining non-freed memory: " << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
            MemoryUsageObserver::getInstance()->purgeInstance();
            logger::purgeInstance();
            errorLogger::purgeInstance();
            return (j == 1 ? 0 : j);
        default:
            break;
    }
    if(molecules->ListOfMolecules.size() == 0){
        mol = new molecule(periode);
        if(mol->cell_size[0] == 0.){
            Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
            for(int i = 0;i < 6;i++){
                Log() << Verbose(1) << "Cell size" << i << ": ";
                cin >> mol->cell_size[i];
            }
        }

        mol->ActiveFlag = true;
        molecules->insert(mol);
    }



    main_menu->perform(molecules, configuration, periode, ConfigFileName);

    if(periode->StorePeriodentafel(configuration->databasepath))
        Log() << Verbose(0) << "Saving of elements.db successful." << endl;

    else
        Log() << Verbose(0) << "Saving of elements.db failed." << endl;

    delete (molecules);
    delete(periode);
  delete(configuration);

  Log() << Verbose(0) <<  "Maximum of allocated memory: "
    << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
  Log() << Verbose(0) <<  "Remaining non-freed memory: "
    << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
  MemoryUsageObserver::purgeInstance();
  logger::purgeInstance();
  errorLogger::purgeInstance();

  return (0);
}

/********************************************** E N D **************************************************/