Changeset 6ac7ee for src/builder.cpp


Ignore:
Timestamp:
Feb 9, 2009, 5:24:10 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
d8b94a
Parents:
124df1
git-author:
Frederik Heber <heber@…> (02/09/09 15:55:37)
git-committer:
Frederik Heber <heber@…> (02/09/09 17:24:10)
Message:

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

Basically, this resulted from a lot of conversions two from spaces to one tab, which is my standard indentation. The mess was caused by eclipse auto-indenting. And in espack2:ConcaveHull was the new stuff, so all from ConcaveHull was replaced in case of doubt.
Additionally, vector had ofstream << operator instead ostream << ...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    • Property mode changed from 100644 to 100755
    r124df1 r6ac7ee  
    1515 * \section about About the Program
    1616 *
    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.
     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.
    2020 *
    21  *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
    22  *  molecular dynamics implementation.
     21 *      A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
     22 *      molecular dynamics implementation.
    2323 *
    2424 * \section install Installation
    2525 *
    26  *  Installation should without problems succeed as follows:
    27  *  -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
    28  *  -# make
    29  *  -# make install
     26 *      Installation should without problems succeed as follows:
     27 *      -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
     28 *      -# make
     29 *      -# make install
    3030 *
    31  *  Further useful commands are
    32  *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
    33  *  -# make doxygen-doc: Creates these html pages out of the documented source
     31 *      Further useful commands are
     32 *      -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
     33 *      -# make doxygen-doc: Creates these html pages out of the documented source
    3434 *
    3535 * \section run Running
    3636 *
    37  *  The program can be executed by running: ./molecuilder
     37 *      The program can be executed by running: ./molecuilder
    3838 *
    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.
     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.
    4242 *
    4343 * \section ref References
    4444 *
    45  *  For the special configuration file format, see the documentation of pcp.
     45 *      For the special configuration file format, see the documentation of pcp.
    4646 *
    4747 */
    … …  
    5050using namespace std;
    5151
     52#include "boundary.hpp"
     53#include "ellipsoid.hpp"
    5254#include "helpers.hpp"
    5355#include "molecules.hpp"
    54 #include "boundary.hpp"
    5556
    5657/********************************************** Submenu routine **************************************/
    … …  
    6263static void AddAtoms(periodentafel *periode, molecule *mol)
    6364{
    64   atom *first, *second, *third, *fourth;
    65   Vector **atoms;
    66   Vector x,y,z,n;  // coordinates for absolute point in cell volume
    67   double a,b,c;
    68   char choice;  // menu choice char
    69   bool valid;
    70 
    71   cout << Verbose(0) << "===========ADD ATOM============================" << endl;
    72   cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
    73   cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
    74   cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
    75   cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
    76   cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
    77   cout << Verbose(0) << "all else - go back" << endl;
    78   cout << Verbose(0) << "===============================================" << endl;
    79   cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
    80   cout << Verbose(0) << "INPUT: ";
    81   cin >> choice;
    82 
    83   switch (choice) {
    84       case 'a': // absolute coordinates of atom
    85         cout << Verbose(0) << "Enter absolute coordinates." << endl;
    86         first = new atom;
    87         first->x.AskPosition(mol->cell_size, false);
    88         first->type = periode->AskElement();  // give type
    89         mol->AddAtom(first);  // add to molecule
    90         break;
    91 
    92       case 'b': // relative coordinates of atom wrt to reference point
    93         first = new atom;
    94         valid = true;
    95         do {
    96           if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
    97           cout << Verbose(0) << "Enter reference coordinates." << endl;
    98           x.AskPosition(mol->cell_size, true);
    99           cout << Verbose(0) << "Enter relative coordinates." << endl;
    100           first->x.AskPosition(mol->cell_size, false);
    101           first->x.AddVector((const Vector *)&x);
    102           cout << Verbose(0) << "\n";
    103         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    104         first->type = periode->AskElement();  // give type
    105         mol->AddAtom(first);  // add to molecule
    106         break;
    107 
    108       case 'c': // relative coordinates of atom wrt to already placed atom
    109         first = new atom;
    110         valid = true;
    111         do {
    112           if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
    113           second = mol->AskAtom("Enter atom number: ");
    114           cout << Verbose(0) << "Enter relative coordinates." << endl;
    115           first->x.AskPosition(mol->cell_size, false);
    116           for (int i=NDIM;i--;) {
    117             first->x.x[i] += second->x.x[i];
    118           }
    119         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    120         first->type = periode->AskElement();  // give type
    121         mol->AddAtom(first);  // add to molecule
    122         break;
    123 
    124       case 'd': // two atoms, two angles and a distance
    125         first = new atom;
    126         valid = true;
    127         do {
    128           if (!valid) {
    129             cout << Verbose(0) << "Resulting coordinates out of cell - ";
    130             first->x.Output((ofstream *)&cout);
    131             cout << Verbose(0) << endl;
    132           }
    133           cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
    134           second = mol->AskAtom("Enter central atom: ");
    135           third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
    136           fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
    137           a = ask_value("Enter distance between central (first) and new atom: ");
    138           b = ask_value("Enter angle between new, first and second atom (degrees): ");
    139           b *= M_PI/180.;
    140           bound(&b, 0., 2.*M_PI);
    141           c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
    142           c *= M_PI/180.;
    143           bound(&c, -M_PI, M_PI);
    144           cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
     65        atom *first, *second, *third, *fourth;
     66        Vector **atoms;
     67        Vector x,y,z,n; // coordinates for absolute point in cell volume
     68        double a,b,c;
     69        char choice;    // menu choice char
     70        bool valid;
     71
     72        cout << Verbose(0) << "===========ADD ATOM============================" << endl;
     73        cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
     74        cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
     75        cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
     76        cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
     77        cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
     78        cout << Verbose(0) << "all else - go back" << endl;
     79        cout << Verbose(0) << "===============================================" << endl;
     80        cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
     81        cout << Verbose(0) << "INPUT: ";
     82        cin >> choice;
     83
     84        switch (choice) {
     85                        case 'a': // absolute coordinates of atom
     86                                cout << Verbose(0) << "Enter absolute coordinates." << endl;
     87                                first = new atom;
     88                                first->x.AskPosition(mol->cell_size, false);
     89                                first->type = periode->AskElement();    // give type
     90                                mol->AddAtom(first);    // add to molecule
     91                                break;
     92
     93                        case 'b': // relative coordinates of atom wrt to reference point
     94                                first = new atom;
     95                                valid = true;
     96                                do {
     97                                        if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
     98                                        cout << Verbose(0) << "Enter reference coordinates." << endl;
     99                                        x.AskPosition(mol->cell_size, true);
     100                                        cout << Verbose(0) << "Enter relative coordinates." << endl;
     101                                        first->x.AskPosition(mol->cell_size, false);
     102                                        first->x.AddVector((const Vector *)&x);
     103                                        cout << Verbose(0) << "\n";
     104                                } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     105                                first->type = periode->AskElement();    // give type
     106                                mol->AddAtom(first);    // add to molecule
     107                                break;
     108
     109                        case 'c': // relative coordinates of atom wrt to already placed atom
     110                                first = new atom;
     111                                valid = true;
     112                                do {
     113                                        if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
     114                                        second = mol->AskAtom("Enter atom number: ");
     115                                        cout << Verbose(0) << "Enter relative coordinates." << endl;
     116                                        first->x.AskPosition(mol->cell_size, false);
     117                                        for (int i=NDIM;i--;) {
     118                                                first->x.x[i] += second->x.x[i];
     119                                        }
     120                                } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     121                                first->type = periode->AskElement();    // give type
     122                                mol->AddAtom(first);    // add to molecule
     123                                break;
     124
     125                        case 'd': // two atoms, two angles and a distance
     126                                first = new atom;
     127                                valid = true;
     128                                do {
     129                                        if (!valid) {
     130                                                cout << Verbose(0) << "Resulting coordinates out of cell - ";
     131                                                first->x.Output((ofstream *)&cout);
     132                                                cout << Verbose(0) << endl;
     133                                        }
     134                                        cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
     135                                        second = mol->AskAtom("Enter central atom: ");
     136                                        third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
     137                                        fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
     138                                        a = ask_value("Enter distance between central (first) and new atom: ");
     139                                        b = ask_value("Enter angle between new, first and second atom (degrees): ");
     140                                        b *= M_PI/180.;
     141                                        bound(&b, 0., 2.*M_PI);
     142                                        c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
     143                                        c *= M_PI/180.;
     144                                        bound(&c, -M_PI, M_PI);
     145                                        cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
    145146/*
    146           second->Output(1,1,(ofstream *)&cout);
    147           third->Output(1,2,(ofstream *)&cout);
    148           fourth->Output(1,3,(ofstream *)&cout);
    149           n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
    150           x.Copyvector(&second->x);
    151           x.SubtractVector(&third->x);
    152           x.Copyvector(&fourth->x);
    153           x.SubtractVector(&third->x);
    154 
    155           if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
    156             cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
    157             continue;
    158           }
    159           cout << Verbose(0) << "resulting relative coordinates: ";
    160           z.Output((ofstream *)&cout);
    161           cout << Verbose(0) << endl;
    162           */
    163           // calc axis vector
    164           x.CopyVector(&second->x);
    165           x.SubtractVector(&third->x);
    166           x.Normalize();
    167           cout << "x: ",
    168           x.Output((ofstream *)&cout);
    169           cout << endl;
    170           z.MakeNormalVector(&second->x,&third->x,&fourth->x);
    171           cout << "z: ",
    172           z.Output((ofstream *)&cout);
    173           cout << endl;
    174           y.MakeNormalVector(&x,&z);
    175           cout << "y: ",
    176           y.Output((ofstream *)&cout);
    177           cout << endl;
    178 
    179           // rotate vector around first angle
    180           first->x.CopyVector(&x);
    181           first->x.RotateVector(&z,b - M_PI);
    182           cout << "Rotated vector: ",
    183           first->x.Output((ofstream *)&cout);
    184           cout << endl;
    185           // remove the projection onto the rotation plane of the second angle
    186           n.CopyVector(&y);
    187           n.Scale(first->x.Projection(&y));
    188           cout << "N1: ",
    189           n.Output((ofstream *)&cout);
    190           cout << endl;
    191           first->x.SubtractVector(&n);
    192           cout << "Subtracted vector: ",
    193           first->x.Output((ofstream *)&cout);
    194           cout << endl;
    195           n.CopyVector(&z);
    196           n.Scale(first->x.Projection(&z));
    197           cout << "N2: ",
    198           n.Output((ofstream *)&cout);
    199           cout << endl;
    200           first->x.SubtractVector(&n);
    201           cout << "2nd subtracted vector: ",
    202           first->x.Output((ofstream *)&cout);
    203           cout << endl;
    204 
    205           // rotate another vector around second angle
    206           n.CopyVector(&y);
    207           n.RotateVector(&x,c - M_PI);
    208           cout << "2nd Rotated vector: ",
    209           n.Output((ofstream *)&cout);
    210           cout << endl;
    211 
    212           // add the two linear independent vectors
    213           first->x.AddVector(&n);
    214           first->x.Normalize();
    215           first->x.Scale(a);
    216           first->x.AddVector(&second->x);
    217 
    218           cout << Verbose(0) << "resulting coordinates: ";
    219           first->x.Output((ofstream *)&cout);
    220           cout << Verbose(0) << endl;
    221         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    222         first->type = periode->AskElement();  // give type
    223         mol->AddAtom(first);  // add to molecule
    224         break;
    225 
    226       case 'e': // least square distance position to a set of atoms
    227         first = new atom;
    228         atoms = new (Vector*[128]);
    229         valid = true;
    230         for(int i=128;i--;)
    231           atoms[i] = NULL;
    232         int i=0, j=0;
    233         cout << Verbose(0) << "Now we need at least three molecules.\n";
    234         do {
    235           cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
    236           cin >> j;
    237           if (j != -1) {
    238             second = mol->FindAtom(j);
    239             atoms[i++] = &(second->x);
    240           }
    241         } while ((j != -1) && (i<128));
    242         if (i >= 2) {
    243           first->x.LSQdistance(atoms, i);
    244 
    245           first->x.Output((ofstream *)&cout);
    246           first->type = periode->AskElement();  // give type
    247           mol->AddAtom(first);  // add to molecule
    248         } else {
    249           delete first;
    250           cout << Verbose(0) << "Please enter at least two vectors!\n";
    251         }
    252         break;
    253   };
     147                                        second->Output(1,1,(ofstream *)&cout);
     148                                        third->Output(1,2,(ofstream *)&cout);
     149                                        fourth->Output(1,3,(ofstream *)&cout);
     150                                        n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
     151                                        x.Copyvector(&second->x);
     152                                        x.SubtractVector(&third->x);
     153                                        x.Copyvector(&fourth->x);
     154                                        x.SubtractVector(&third->x);
     155
     156                                        if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
     157                                                cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
     158                                                continue;
     159                                        }
     160                                        cout << Verbose(0) << "resulting relative coordinates: ";
     161                                        z.Output((ofstream *)&cout);
     162                                        cout << Verbose(0) << endl;
     163                                        */
     164                                        // calc axis vector
     165                                        x.CopyVector(&second->x);
     166                                        x.SubtractVector(&third->x);
     167                                        x.Normalize();
     168                                        cout << "x: ",
     169                                        x.Output((ofstream *)&cout);
     170                                        cout << endl;
     171                                        z.MakeNormalVector(&second->x,&third->x,&fourth->x);
     172                                        cout << "z: ",
     173                                        z.Output((ofstream *)&cout);
     174                                        cout << endl;
     175                                        y.MakeNormalVector(&x,&z);
     176                                        cout << "y: ",
     177                                        y.Output((ofstream *)&cout);
     178                                        cout << endl;
     179
     180                                        // rotate vector around first angle
     181                                        first->x.CopyVector(&x);
     182                                        first->x.RotateVector(&z,b - M_PI);
     183                                        cout << "Rotated vector: ",
     184                                        first->x.Output((ofstream *)&cout);
     185                                        cout << endl;
     186                                        // remove the projection onto the rotation plane of the second angle
     187                                        n.CopyVector(&y);
     188                                        n.Scale(first->x.Projection(&y));
     189                                        cout << "N1: ",
     190                                        n.Output((ofstream *)&cout);
     191                                        cout << endl;
     192                                        first->x.SubtractVector(&n);
     193                                        cout << "Subtracted vector: ",
     194                                        first->x.Output((ofstream *)&cout);
     195                                        cout << endl;
     196                                        n.CopyVector(&z);
     197                                        n.Scale(first->x.Projection(&z));
     198                                        cout << "N2: ",
     199                                        n.Output((ofstream *)&cout);
     200                                        cout << endl;
     201                                        first->x.SubtractVector(&n);
     202                                        cout << "2nd subtracted vector: ",
     203                                        first->x.Output((ofstream *)&cout);
     204                                        cout << endl;
     205
     206                                        // rotate another vector around second angle
     207                                        n.CopyVector(&y);
     208                                        n.RotateVector(&x,c - M_PI);
     209                                        cout << "2nd Rotated vector: ",
     210                                        n.Output((ofstream *)&cout);
     211                                        cout << endl;
     212
     213                                        // add the two linear independent vectors
     214                                        first->x.AddVector(&n);
     215                                        first->x.Normalize();
     216                                        first->x.Scale(a);
     217                                        first->x.AddVector(&second->x);
     218
     219                                        cout << Verbose(0) << "resulting coordinates: ";
     220                                        first->x.Output((ofstream *)&cout);
     221                                        cout << Verbose(0) << endl;
     222                                } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     223                                first->type = periode->AskElement();    // give type
     224                                mol->AddAtom(first);    // add to molecule
     225                                break;
     226
     227                        case 'e': // least square distance position to a set of atoms
     228                                first = new atom;
     229                                atoms = new (Vector*[128]);
     230                                valid = true;
     231                                for(int i=128;i--;)
     232                                        atoms[i] = NULL;
     233                                int i=0, j=0;
     234                                cout << Verbose(0) << "Now we need at least three molecules.\n";
     235                                do {
     236                                        cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
     237                                        cin >> j;
     238                                        if (j != -1) {
     239                                                second = mol->FindAtom(j);
     240                                                atoms[i++] = &(second->x);
     241                                        }
     242                                } while ((j != -1) && (i<128));
     243                                if (i >= 2) {
     244                                        first->x.LSQdistance(atoms, i);
     245
     246                                        first->x.Output((ofstream *)&cout);
     247                                        first->type = periode->AskElement();    // give type
     248                                        mol->AddAtom(first);    // add to molecule
     249                                } else {
     250                                        delete first;
     251                                        cout << Verbose(0) << "Please enter at least two vectors!\n";
     252                                }
     253                                break;
     254        };
    254255};
    255256
    … …  
    259260static void CenterAtoms(molecule *mol)
    260261{
    261   Vector x, y;
    262   char choice;  // menu choice char
    263 
    264   cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
    265   cout << Verbose(0) << " a - on origin" << endl;
    266   cout << Verbose(0) << " b - on center of gravity" << endl;
    267   cout << Verbose(0) << " c - within box with additional boundary" << endl;
    268   cout << Verbose(0) << " d - within given simulation box" << endl;
    269   cout << Verbose(0) << "all else - go back" << endl;
    270   cout << Verbose(0) << "===============================================" << endl;
    271   cout << Verbose(0) << "INPUT: ";
    272   cin >> choice;
    273 
    274   switch (choice) {
    275     default:
    276       cout << Verbose(0) << "Not a valid choice." << endl;
    277       break;
    278     case 'a':
    279       cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
    280       mol->CenterOrigin((ofstream *)&cout, &x);
    281       break;
    282     case 'b':
    283       cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
    284       mol->CenterGravity((ofstream *)&cout, &x);
    285       break;
    286     case 'c':
    287       cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
    288       for (int i=0;i<NDIM;i++) {
    289         cout << Verbose(0) << "Enter axis " << i << " boundary: ";
    290         cin >> y.x[i];
    291       }
    292       mol->CenterEdge((ofstream *)&cout, &x);  // make every coordinate positive
    293       mol->Translate(&y); // translate by boundary
    294       mol->SetBoxDimension(&(x+y*2));  // update Box of atoms by boundary
    295       break;
    296     case 'd':
    297       cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
    298       for (int i=0;i<NDIM;i++) {
    299         cout << Verbose(0) << "Enter axis " << i << " boundary: ";
    300         cin >> x.x[i];
    301       }
    302       // center
    303       mol->CenterInBox((ofstream *)&cout, &x);
    304       // update Box of atoms by boundary
    305       mol->SetBoxDimension(&x);
    306       break;
    307   }
     262        Vector x, y, helper;
     263        char choice;    // menu choice char
     264
     265        cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
     266        cout << Verbose(0) << " a - on origin" << endl;
     267        cout << Verbose(0) << " b - on center of gravity" << endl;
     268        cout << Verbose(0) << " c - within box with additional boundary" << endl;
     269        cout << Verbose(0) << " d - within given simulation box" << endl;
     270        cout << Verbose(0) << "all else - go back" << endl;
     271        cout << Verbose(0) << "===============================================" << endl;
     272        cout << Verbose(0) << "INPUT: ";
     273        cin >> choice;
     274
     275        switch (choice) {
     276                default:
     277                        cout << Verbose(0) << "Not a valid choice." << endl;
     278                        break;
     279                case 'a':
     280                        cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
     281                        mol->CenterOrigin((ofstream *)&cout, &x);
     282                        break;
     283                case 'b':
     284                        cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
     285                        mol->CenterGravity((ofstream *)&cout, &x);
     286                        break;
     287                case 'c':
     288                        cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
     289                        for (int i=0;i<NDIM;i++) {
     290                                cout << Verbose(0) << "Enter axis " << i << " boundary: ";
     291                                cin >> y.x[i];
     292                        }
     293                        mol->CenterEdge((ofstream *)&cout, &x); // make every coordinate positive
     294                        mol->Translate(&y); // translate by boundary
     295                        helper.CopyVector(&y);
     296                        helper.Scale(2.);
     297                        helper.AddVector(&x);
     298                        mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
     299                        break;
     300                case 'd':
     301                        cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     302                        for (int i=0;i<NDIM;i++) {
     303                                cout << Verbose(0) << "Enter axis " << i << " boundary: ";
     304                                cin >> x.x[i];
     305                        }
     306                        // center
     307                        mol->CenterInBox((ofstream *)&cout, &x);
     308                        // update Box of atoms by boundary
     309                        mol->SetBoxDimension(&x);
     310                        break;
     311        }
    308312};
    309313
    … …  
    314318static void AlignAtoms(periodentafel *periode, molecule *mol)
    315319{
    316   atom *first, *second, *third;
    317   Vector x,n;
    318   char choice;  // menu choice char
    319 
    320   cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
    321   cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
    322   cout << Verbose(0) << " b - state alignment vector" << endl;
    323   cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
    324   cout << Verbose(0) << " d - align automatically by least square fit" << endl;
    325   cout << Verbose(0) << "all else - go back" << endl;
    326   cout << Verbose(0) << "===============================================" << endl;
    327   cout << Verbose(0) << "INPUT: ";
    328   cin >> choice;
    329 
    330   switch (choice) {
    331     default:
    332     case 'a': // three atoms defining mirror plane
    333       first = mol->AskAtom("Enter first atom: ");
    334       second = mol->AskAtom("Enter second atom: ");
    335       third = mol->AskAtom("Enter third atom: ");
    336 
    337       n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
    338       break;
    339     case 'b': // normal vector of mirror plane
    340       cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    341       n.AskPosition(mol->cell_size,0);
    342       n.Normalize();
    343       break;
    344     case 'c': // three atoms defining mirror plane
    345       first = mol->AskAtom("Enter first atom: ");
    346       second = mol->AskAtom("Enter second atom: ");
    347 
    348       n.CopyVector((const Vector *)&first->x);
    349       n.SubtractVector((const Vector *)&second->x);
    350       n.Normalize();
    351       break;
    352     case 'd':
    353       char shorthand[4];
    354       Vector a;
    355       struct lsq_params param;
    356       do {
    357         fprintf(stdout, "Enter the element of atoms to be chosen: ");
    358         fscanf(stdin, "%3s", shorthand);
    359       } while ((param.type = periode->FindElement(shorthand)) == NULL);
    360       cout << Verbose(0) << "Element is " << param.type->name << endl;
    361       mol->GetAlignvector(&param);
    362       for (int i=NDIM;i--;) {
    363         x.x[i] = gsl_vector_get(param.x,i);
    364         n.x[i] = gsl_vector_get(param.x,i+NDIM);
    365       }
    366       gsl_vector_free(param.x);
    367       cout << Verbose(0) << "Offset vector: ";
    368       x.Output((ofstream *)&cout);
    369       cout << Verbose(0) << endl;
    370       n.Normalize();
    371       break;
    372   };
    373   cout << Verbose(0) << "Alignment vector: ";
    374   n.Output((ofstream *)&cout);
    375   cout << Verbose(0) << endl;
    376   mol->Align(&n);
     320        atom *first, *second, *third;
     321        Vector x,n;
     322        char choice;    // menu choice char
     323
     324        cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
     325        cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
     326        cout << Verbose(0) << " b - state alignment vector" << endl;
     327        cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
     328        cout << Verbose(0) << " d - align automatically by least square fit" << endl;
     329        cout << Verbose(0) << "all else - go back" << endl;
     330        cout << Verbose(0) << "===============================================" << endl;
     331        cout << Verbose(0) << "INPUT: ";
     332        cin >> choice;
     333
     334        switch (choice) {
     335                default:
     336                case 'a': // three atoms defining mirror plane
     337                        first = mol->AskAtom("Enter first atom: ");
     338                        second = mol->AskAtom("Enter second atom: ");
     339                        third = mol->AskAtom("Enter third atom: ");
     340
     341                        n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
     342                        break;
     343                case 'b': // normal vector of mirror plane
     344                        cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
     345                        n.AskPosition(mol->cell_size,0);
     346                        n.Normalize();
     347                        break;
     348                case 'c': // three atoms defining mirror plane
     349                        first = mol->AskAtom("Enter first atom: ");
     350                        second = mol->AskAtom("Enter second atom: ");
     351
     352                        n.CopyVector((const Vector *)&first->x);
     353                        n.SubtractVector((const Vector *)&second->x);
     354                        n.Normalize();
     355                        break;
     356                case 'd':
     357                        char shorthand[4];
     358                        Vector a;
     359                        struct lsq_params param;
     360                        do {
     361                                fprintf(stdout, "Enter the element of atoms to be chosen: ");
     362                                fscanf(stdin, "%3s", shorthand);
     363                        } while ((param.type = periode->FindElement(shorthand)) == NULL);
     364                        cout << Verbose(0) << "Element is " << param.type->name << endl;
     365                        mol->GetAlignvector(&param);
     366                        for (int i=NDIM;i--;) {
     367                                x.x[i] = gsl_vector_get(param.x,i);
     368                                n.x[i] = gsl_vector_get(param.x,i+NDIM);
     369                        }
     370                        gsl_vector_free(param.x);
     371                        cout << Verbose(0) << "Offset vector: ";
     372                        x.Output((ofstream *)&cout);
     373                        cout << Verbose(0) << endl;
     374                        n.Normalize();
     375                        break;
     376        };
     377        cout << Verbose(0) << "Alignment vector: ";
     378        n.Output((ofstream *)&cout);
     379        cout << Verbose(0) << endl;
     380        mol->Align(&n);
    377381};
    378382
    … …  
    382386static void MirrorAtoms(molecule *mol)
    383387{
    384   atom *first, *second, *third;
    385   Vector n;
    386   char choice;  // menu choice char
    387 
    388   cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
    389   cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
    390   cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
    391   cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
    392   cout << Verbose(0) << "all else - go back" << endl;
    393   cout << Verbose(0) << "===============================================" << endl;
    394   cout << Verbose(0) << "INPUT: ";
    395   cin >> choice;
    396 
    397   switch (choice) {
    398     default:
    399     case 'a': // three atoms defining mirror plane
    400       first = mol->AskAtom("Enter first atom: ");
    401       second = mol->AskAtom("Enter second atom: ");
    402       third = mol->AskAtom("Enter third atom: ");
    403 
    404       n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
    405       break;
    406     case 'b': // normal vector of mirror plane
    407       cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    408       n.AskPosition(mol->cell_size,0);
    409       n.Normalize();
    410       break;
    411     case 'c': // three atoms defining mirror plane
    412       first = mol->AskAtom("Enter first atom: ");
    413       second = mol->AskAtom("Enter second atom: ");
    414 
    415       n.CopyVector((const Vector *)&first->x);
    416       n.SubtractVector((const Vector *)&second->x);
    417       n.Normalize();
    418       break;
    419   };
    420   cout << Verbose(0) << "Normal vector: ";
    421   n.Output((ofstream *)&cout);
    422   cout << Verbose(0) << endl;
    423   mol->Mirror((const Vector *)&n);
     388        atom *first, *second, *third;
     389        Vector n;
     390        char choice;    // menu choice char
     391
     392        cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
     393        cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
     394        cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
     395        cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
     396        cout << Verbose(0) << "all else - go back" << endl;
     397        cout << Verbose(0) << "===============================================" << endl;
     398        cout << Verbose(0) << "INPUT: ";
     399        cin >> choice;
     400
     401        switch (choice) {
     402                default:
     403                case 'a': // three atoms defining mirror plane
     404                        first = mol->AskAtom("Enter first atom: ");
     405                        second = mol->AskAtom("Enter second atom: ");
     406                        third = mol->AskAtom("Enter third atom: ");
     407
     408                        n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
     409                        break;
     410                case 'b': // normal vector of mirror plane
     411                        cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
     412                        n.AskPosition(mol->cell_size,0);
     413                        n.Normalize();
     414                        break;
     415                case 'c': // three atoms defining mirror plane
     416                        first = mol->AskAtom("Enter first atom: ");
     417                        second = mol->AskAtom("Enter second atom: ");
     418
     419                        n.CopyVector((const Vector *)&first->x);
     420                        n.SubtractVector((const Vector *)&second->x);
     421                        n.Normalize();
     422                        break;
     423        };
     424        cout << Verbose(0) << "Normal vector: ";
     425        n.Output((ofstream *)&cout);
     426        cout << Verbose(0) << endl;
     427        mol->Mirror((const Vector *)&n);
    424428};
    425429
    … …  
    429433static void RemoveAtoms(molecule *mol)
    430434{
    431   atom *first, *second;
    432   int axis;
    433   double tmp1, tmp2;
    434   char choice;  // menu choice char
    435 
    436   cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
    437   cout << Verbose(0) << " a - state atom for removal by number" << endl;
    438   cout << Verbose(0) << " b - keep only in radius around atom" << endl;
    439   cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
    440   cout << Verbose(0) << "all else - go back" << endl;
    441   cout << Verbose(0) << "===============================================" << endl;
    442   cout << Verbose(0) << "INPUT: ";
    443   cin >> choice;
    444 
    445   switch (choice) {
    446     default:
    447     case 'a':
    448       if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
    449         cout << Verbose(1) << "Atom removed." << endl;
    450       else
    451         cout << Verbose(1) << "Atom not found." << endl;
    452       break;
    453     case 'b':
    454       second = mol->AskAtom("Enter number of atom as reference point: ");
    455       cout << Verbose(0) << "Enter radius: ";
    456       cin >> tmp1;
    457       first = mol->start;
    458       while(first->next != mol->end) {
    459         first = first->next;
    460         if (first->x.Distance((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
    461           mol->RemoveAtom(first);
    462       }
    463       break;
    464     case 'c':
    465       cout << Verbose(0) << "Which axis is it: ";
    466       cin >> axis;
    467       cout << Verbose(0) << "Left inward boundary: ";
    468       cin >> tmp1;
    469       cout << Verbose(0) << "Right inward boundary: ";
    470       cin >> tmp2;
    471       first = mol->start;
    472       while(first->next != mol->end) {
    473         first = first->next;
    474         if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...
    475           mol->RemoveAtom(first);
    476       }
    477       break;
    478   };
    479   //mol->Output((ofstream *)&cout);
    480   choice = 'r';
     435        atom *first, *second;
     436        int axis;
     437        double tmp1, tmp2;
     438        char choice;    // menu choice char
     439
     440        cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
     441        cout << Verbose(0) << " a - state atom for removal by number" << endl;
     442        cout << Verbose(0) << " b - keep only in radius around atom" << endl;
     443        cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
     444        cout << Verbose(0) << "all else - go back" << endl;
     445        cout << Verbose(0) << "===============================================" << endl;
     446        cout << Verbose(0) << "INPUT: ";
     447        cin >> choice;
     448
     449        switch (choice) {
     450                default:
     451                case 'a':
     452                        if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
     453                                cout << Verbose(1) << "Atom removed." << endl;
     454                        else
     455                                cout << Verbose(1) << "Atom not found." << endl;
     456                        break;
     457                case 'b':
     458                        second = mol->AskAtom("Enter number of atom as reference point: ");
     459                        cout << Verbose(0) << "Enter radius: ";
     460                        cin >> tmp1;
     461                        first = mol->start;
     462                        while(first->next != mol->end) {
     463                                first = first->next;
     464                                if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
     465                                        mol->RemoveAtom(first);
     466                        }
     467                        break;
     468                case 'c':
     469                        cout << Verbose(0) << "Which axis is it: ";
     470                        cin >> axis;
     471                        cout << Verbose(0) << "Left inward boundary: ";
     472                        cin >> tmp1;
     473                        cout << Verbose(0) << "Right inward boundary: ";
     474                        cin >> tmp2;
     475                        first = mol->start;
     476                        while(first->next != mol->end) {
     477                                first = first->next;
     478                                if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...
     479                                        mol->RemoveAtom(first);
     480                        }
     481                        break;
     482        };
     483        //mol->Output((ofstream *)&cout);
     484        choice = 'r';
    481485};
    482486
    … …  
    487491static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
    488492{
    489   atom *first, *second, *third;
    490   Vector x,y;
    491   double min[256], tmp1, tmp2, tmp3;
    492   int Z;
    493   char choice;  // menu choice char
    494 
    495   cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
    496   cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
    497   cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
    498   cout << Verbose(0) << " c - calculate bond angle" << endl;
    499   cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
    500   cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
    501   cout << Verbose(0) << " f - calculate temperature from current velocity" << endl;
    502   cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
    503   cout << Verbose(0) << "all else - go back" << endl;
    504   cout << Verbose(0) << "===============================================" << endl;
    505   cout << Verbose(0) << "INPUT: ";
    506   cin >> choice;
    507 
    508   switch(choice) {
    509     default:
    510       cout << Verbose(1) << "Not a valid choice." << endl;
    511       break;
    512     case 'a':
    513       first = mol->AskAtom("Enter first atom: ");
    514       for (int i=MAX_ELEMENTS;i--;)
    515         min[i] = 0.;
    516 
    517       second = mol->start;
    518       while ((second->next != mol->end)) {
    519         second = second->next; // advance
    520         Z = second->type->Z;
    521         tmp1 = 0.;
    522         if (first != second) {
    523           x.CopyVector((const Vector *)&first->x);
    524           x.SubtractVector((const Vector *)&second->x);
    525           tmp1 = x.Norm();
    526         }
    527         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
    528         //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
    529       }
    530       for (int i=MAX_ELEMENTS;i--;)
    531         if (min[i] != 0.) cout << 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;
    532       break;
    533 
    534     case 'b':
    535       first = mol->AskAtom("Enter first atom: ");
    536       second = mol->AskAtom("Enter second atom: ");
    537       for (int i=NDIM;i--;)
    538         min[i] = 0.;
    539       x.CopyVector((const Vector *)&first->x);
    540       x.SubtractVector((const Vector *)&second->x);
    541       tmp1 = x.Norm();
    542       cout << Verbose(1) << "Distance vector is ";
    543       x.Output((ofstream *)&cout);
    544       cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
    545       break;
    546 
    547     case 'c':
    548       cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
    549       first = mol->AskAtom("Enter first atom: ");
    550       second = mol->AskAtom("Enter central atom: ");
    551       third  = mol->AskAtom("Enter last atom: ");
    552       tmp1 = tmp2 = tmp3 = 0.;
    553       x.CopyVector((const Vector *)&first->x);
    554       x.SubtractVector((const Vector *)&second->x);
    555       y.CopyVector((const Vector *)&third->x);
    556       y.SubtractVector((const Vector *)&second->x);
    557       cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
    558       cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
    559       break;
    560     case 'd':
    561         cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
    562         cout << Verbose(0) << "Shall we rotate? [0/1]: ";
    563         cin >> Z;
    564         if ((Z >=0) && (Z <=1))
    565           mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
    566         else
    567           mol->PrincipalAxisSystem((ofstream *)&cout, false);
    568         break;
    569     case 'e':
    570         cout << Verbose(0) << "Evaluating volume of the convex envelope.";
    571         VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
    572         break;
    573     case 'f':
    574       mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
    575       break;
    576     case 'g':
    577       {
    578         char filename[255];
    579         cout << "Please enter filename: " << endl;
    580         cin >> filename;
    581         cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
    582         ofstream *output = new ofstream(filename, ios::trunc);
    583         if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
    584           cout << Verbose(2) << "File could not be written." << endl;
    585         else
    586           cout << Verbose(2) << "File stored." << endl;
    587         output->close();
    588         delete(output);
    589       }
    590       break;
    591   }
     493        atom *first, *second, *third;
     494        Vector x,y;
     495        double min[256], tmp1, tmp2, tmp3;
     496        int Z;
     497        char choice;    // menu choice char
     498
     499        cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
     500        cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
     501        cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
     502        cout << Verbose(0) << " c - calculate bond angle" << endl;
     503        cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
     504        cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
     505        cout << Verbose(0) << " f - calculate temperature from current velocity" << endl;
     506        cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
     507        cout << Verbose(0) << "all else - go back" << endl;
     508        cout << Verbose(0) << "===============================================" << endl;
     509        cout << Verbose(0) << "INPUT: ";
     510        cin >> choice;
     511
     512        switch(choice) {
     513                default:
     514                        cout << Verbose(1) << "Not a valid choice." << endl;
     515                        break;
     516                case 'a':
     517                        first = mol->AskAtom("Enter first atom: ");
     518                        for (int i=MAX_ELEMENTS;i--;)
     519                                min[i] = 0.;
     520
     521                        second = mol->start;
     522                        while ((second->next != mol->end)) {
     523                                second = second->next; // advance
     524                                Z = second->type->Z;
     525                                tmp1 = 0.;
     526                                if (first != second) {
     527                                        x.CopyVector((const Vector *)&first->x);
     528                                        x.SubtractVector((const Vector *)&second->x);
     529                                        tmp1 = x.Norm();
     530                                }
     531                                if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
     532                                //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
     533                        }
     534                        for (int i=MAX_ELEMENTS;i--;)
     535                                if (min[i] != 0.) cout << 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;
     536                        break;
     537
     538                case 'b':
     539                        first = mol->AskAtom("Enter first atom: ");
     540                        second = mol->AskAtom("Enter second atom: ");
     541                        for (int i=NDIM;i--;)
     542                                min[i] = 0.;
     543                        x.CopyVector((const Vector *)&first->x);
     544                        x.SubtractVector((const Vector *)&second->x);
     545                        tmp1 = x.Norm();
     546                        cout << Verbose(1) << "Distance vector is ";
     547                        x.Output((ofstream *)&cout);
     548                        cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
     549                        break;
     550
     551                case 'c':
     552                        cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
     553                        first = mol->AskAtom("Enter first atom: ");
     554                        second = mol->AskAtom("Enter central atom: ");
     555                        third   = mol->AskAtom("Enter last atom: ");
     556                        tmp1 = tmp2 = tmp3 = 0.;
     557                        x.CopyVector((const Vector *)&first->x);
     558                        x.SubtractVector((const Vector *)&second->x);
     559                        y.CopyVector((const Vector *)&third->x);
     560                        y.SubtractVector((const Vector *)&second->x);
     561                        cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     562                        cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
     563                        break;
     564                case 'd':
     565                        cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
     566                        cout << Verbose(0) << "Shall we rotate? [0/1]: ";
     567                        cin >> Z;
     568                        if ((Z >=0) && (Z <=1))
     569                                mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
     570                        else
     571                                mol->PrincipalAxisSystem((ofstream *)&cout, false);
     572                        break;
     573                case 'e':
     574                        cout << Verbose(0) << "Evaluating volume of the convex envelope.";
     575                        VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
     576                        break;
     577                case 'f':
     578                        mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
     579                        break;
     580                case 'g':
     581                        {
     582                                char filename[255];
     583                                cout << "Please enter filename: " << endl;
     584                                cin >> filename;
     585                                cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
     586                                ofstream *output = new ofstream(filename, ios::trunc);
     587                                if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
     588                                        cout << Verbose(2) << "File could not be written." << endl;
     589                                else
     590                                        cout << Verbose(2) << "File stored." << endl;
     591                                output->close();
     592                                delete(output);
     593                        }
     594                        break;
     595        }
    592596};
    593597
    … …  
    598602static void FragmentAtoms(molecule *mol, config *configuration)
    599603{
    600   int Order1;
    601   clock_t start, end;
    602 
    603   cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
    604   cout << Verbose(0) << "What's the desired bond order: ";
    605   cin >> Order1;
    606   if (mol->first->next != mol->last) {  // there are bonds
    607     start = clock();
    608     mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);
    609     end = clock();
    610     cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    611   } else
    612     cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
     604        int Order1;
     605        clock_t start, end;
     606
     607        cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
     608        cout << Verbose(0) << "What's the desired bond order: ";
     609        cin >> Order1;
     610        if (mol->first->next != mol->last) {    // there are bonds
     611                start = clock();
     612                mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);
     613                end = clock();
     614                cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     615        } else
     616                cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
    613617};
    614618
    … …  
    619623static void testroutine(molecule *mol)
    620624{
    621   // the current test routine checks the functionality of the KeySet&Graph concept:
    622   // We want to have a multiindex (the KeySet) describing a unique subgraph
    623   atom *Walker = mol->start;
    624   int i, comp, counter=0;
    625 
    626   // generate some KeySets
    627   cout << "Generating KeySets." << endl;
    628   KeySet TestSets[mol->AtomCount+1];
    629   i=1;
    630   while (Walker->next != mol->end) {
    631     Walker = Walker->next;
    632     for (int j=0;j<i;j++) {
    633       TestSets[j].insert(Walker->nr);
    634     }
    635     i++;
    636   }
    637   cout << "Testing insertion of already present item in KeySets." << endl;
    638   KeySetTestPair test;
    639   test = TestSets[mol->AtomCount-1].insert(Walker->nr);
    640   if (test.second) {
    641     cout << Verbose(1) << "Insertion worked?!" << endl;
    642   } else {
    643     cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
    644   }
    645   TestSets[mol->AtomCount].insert(mol->end->previous->nr);
    646   TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
    647 
    648   // constructing Graph structure
    649   cout << "Generating Subgraph class." << endl;
    650   Graph Subgraphs;
    651 
    652   // insert KeySets into Subgraphs
    653   cout << "Inserting KeySets into Subgraph class." << endl;
    654   for (int j=0;j<mol->AtomCount;j++) {
    655     Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
    656   }
    657   cout << "Testing insertion of already present item in Subgraph." << endl;
    658   GraphTestPair test2;
    659   test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
    660   if (test2.second) {
    661     cout << Verbose(1) << "Insertion worked?!" << endl;
    662   } else {
    663     cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
    664   }
    665 
    666   // show graphs
    667   cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
    668   Graph::iterator A = Subgraphs.begin();
    669   while (A !=  Subgraphs.end()) {
    670     cout << (*A).second.first << ": ";
    671     KeySet::iterator key = (*A).first.begin();
    672     comp = -1;
    673     while (key != (*A).first.end()) {
    674       if ((*key) > comp)
    675         cout << (*key) << " ";
    676       else
    677         cout << (*key) << "! ";
    678       comp = (*key);
    679       key++;
    680     }
    681     cout << endl;
    682     A++;
    683   }
     625        // the current test routine checks the functionality of the KeySet&Graph concept:
     626        // We want to have a multiindex (the KeySet) describing a unique subgraph
     627        atom *Walker = mol->start;
     628        int i, comp, counter=0;
     629
     630        // generate some KeySets
     631        cout << "Generating KeySets." << endl;
     632        KeySet TestSets[mol->AtomCount+1];
     633        i=1;
     634        while (Walker->next != mol->end) {
     635                Walker = Walker->next;
     636                for (int j=0;j<i;j++) {
     637                        TestSets[j].insert(Walker->nr);
     638                }
     639                i++;
     640        }
     641        cout << "Testing insertion of already present item in KeySets." << endl;
     642        KeySetTestPair test;
     643        test = TestSets[mol->AtomCount-1].insert(Walker->nr);
     644        if (test.second) {
     645                cout << Verbose(1) << "Insertion worked?!" << endl;
     646        } else {
     647                cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
     648        }
     649        TestSets[mol->AtomCount].insert(mol->end->previous->nr);
     650        TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
     651
     652        // constructing Graph structure
     653        cout << "Generating Subgraph class." << endl;
     654        Graph Subgraphs;
     655
     656        // insert KeySets into Subgraphs
     657        cout << "Inserting KeySets into Subgraph class." << endl;
     658        for (int j=0;j<mol->AtomCount;j++) {
     659                Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
     660        }
     661        cout << "Testing insertion of already present item in Subgraph." << endl;
     662        GraphTestPair test2;
     663        test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
     664        if (test2.second) {
     665                cout << Verbose(1) << "Insertion worked?!" << endl;
     666        } else {
     667                cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
     668        }
     669
     670        // show graphs
     671        cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
     672        Graph::iterator A = Subgraphs.begin();
     673        while (A !=     Subgraphs.end()) {
     674                cout << (*A).second.first << ": ";
     675                KeySet::iterator key = (*A).first.begin();
     676                comp = -1;
     677                while (key != (*A).first.end()) {
     678                        if ((*key) > comp)
     679                                cout << (*key) << " ";
     680                        else
     681                                cout << (*key) << "! ";
     682                        comp = (*key);
     683                        key++;
     684                }
     685                cout << endl;
     686                A++;
     687        }
    684688};
    685689
    … …  
    692696static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, molecule *mol)
    693697{
    694   char filename[MAXSTRINGSIZE];
    695   ofstream output;
    696 
    697   cout << Verbose(0) << "Storing configuration ... " << endl;
    698   // get correct valence orbitals
    699   mol->CalculateOrbitals(*configuration);
    700   configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
    701   strcpy(filename, ConfigFileName);
    702   if (ConfigFileName != NULL) { // test the file name
    703     output.open(ConfigFileName, ios::trunc);
    704   } else if (strlen(configuration->configname) != 0) {
    705     strcpy(filename, configuration->configname);
    706     output.open(configuration->configname, ios::trunc);
    707     } else {
    708       strcpy(filename, DEFAULTCONFIG);
    709       output.open(DEFAULTCONFIG, ios::trunc);
    710     }
    711   output.close();
    712   output.clear();
    713   cout << Verbose(0) << "Saving of config file ";
    714   if (configuration->Save(filename, periode, mol))
    715     cout << "successful." << endl;
    716   else
    717     cout << "failed." << endl;
    718 
    719   // and save to xyz file
    720   if (ConfigFileName != NULL) {
    721     strcpy(filename, ConfigFileName);
    722     strcat(filename, ".xyz");
    723     output.open(filename, ios::trunc);
    724   }
    725   if (output == NULL) {
    726     strcpy(filename,"main_pcp_linux");
    727     strcat(filename, ".xyz");
    728     output.open(filename, ios::trunc);
    729   }
    730   cout << Verbose(0) << "Saving of XYZ file ";
    731   if (mol->MDSteps <= 1) {
    732     if (mol->OutputXYZ(&output))
    733       cout << "successful." << endl;
    734     else
    735       cout << "failed." << endl;
    736   } else {
    737     if (mol->OutputTrajectoriesXYZ(&output))
    738       cout << "successful." << endl;
    739     else
    740       cout << "failed." << endl;
    741   }
    742   output.close();
    743   output.clear();
    744 
    745   // and save as MPQC configuration
    746   if (ConfigFileName != NULL)
    747     strcpy(filename, ConfigFileName);
    748   if (output == NULL)
    749     strcpy(filename,"main_pcp_linux");
    750   cout << Verbose(0) << "Saving as mpqc input ";
    751   if (configuration->SaveMPQC(filename, mol))
    752     cout << "done." << endl;
    753   else
    754     cout << "failed." << endl;
    755 
    756   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    757     cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
    758   }
     698        char filename[MAXSTRINGSIZE];
     699        ofstream output;
     700
     701        cout << Verbose(0) << "Storing configuration ... " << endl;
     702        // get correct valence orbitals
     703        mol->CalculateOrbitals(*configuration);
     704        configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
     705        strcpy(filename, ConfigFileName);
     706        if (ConfigFileName != NULL) { // test the file name
     707                output.open(ConfigFileName, ios::trunc);
     708        } else if (strlen(configuration->configname) != 0) {
     709                strcpy(filename, configuration->configname);
     710                output.open(configuration->configname, ios::trunc);
     711                } else {
     712                        strcpy(filename, DEFAULTCONFIG);
     713                        output.open(DEFAULTCONFIG, ios::trunc);
     714                }
     715        output.close();
     716        output.clear();
     717        cout << Verbose(0) << "Saving of config file ";
     718        if (configuration->Save(filename, periode, mol))
     719                cout << "successful." << endl;
     720        else
     721                cout << "failed." << endl;
     722
     723        // and save to xyz file
     724        if (ConfigFileName != NULL) {
     725                strcpy(filename, ConfigFileName);
     726                strcat(filename, ".xyz");
     727                output.open(filename, ios::trunc);
     728        }
     729        if (output == NULL) {
     730                strcpy(filename,"main_pcp_linux");
     731                strcat(filename, ".xyz");
     732                output.open(filename, ios::trunc);
     733        }
     734        cout << Verbose(0) << "Saving of XYZ file ";
     735        if (mol->MDSteps <= 1) {
     736                if (mol->OutputXYZ(&output))
     737                        cout << "successful." << endl;
     738                else
     739                        cout << "failed." << endl;
     740        } else {
     741                if (mol->OutputTrajectoriesXYZ(&output))
     742                        cout << "successful." << endl;
     743                else
     744                        cout << "failed." << endl;
     745        }
     746        output.close();
     747        output.clear();
     748
     749        // and save as MPQC configuration
     750        if (ConfigFileName != NULL)
     751                strcpy(filename, ConfigFileName);
     752        if (output == NULL)
     753                strcpy(filename,"main_pcp_linux");
     754        cout << Verbose(0) << "Saving as mpqc input ";
     755        if (configuration->SaveMPQC(filename, mol))
     756                cout << "done." << endl;
     757        else
     758                cout << "failed." << endl;
     759
     760        if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
     761                cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
     762        }
    759763};
    760764
    … …  
    771775static int ParseCommandLineOptions(int argc, char **argv, molecule *&mol, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
    772776{
    773   Vector x,y,z,n;  // coordinates for absolute point in cell volume
    774   double *factor; // unit factor if desired
    775   ifstream test;
    776   ofstream output;
    777   string line;
    778   atom *first;
    779   bool SaveFlag = false;
    780   int ExitFlag = 0;
    781   int j;
    782   double volume = 0.;
    783   enum ConfigStatus config_present = absent;
    784   clock_t start,end;
    785   int argptr;
    786   PathToDatabases = LocalPath;
    787 
    788   if (argc > 1) { // config file specified as option
    789     // 1. : Parse options that just set variables or print help
    790     argptr = 1;
    791     do {
    792       if (argv[argptr][0] == '-') {
    793         cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
    794         argptr++;
    795         switch(argv[argptr-1][1]) {
    796           case 'h':
    797           case 'H':
    798           case '?':
    799             cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
    800             cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
    801             cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
    802             cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
    803             cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
    804             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
    805             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    806             cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    807             cout << "\t-O\tCenter atoms in origin." << endl;
    808             cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    809             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    810             cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
    811             cout << "\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;
    812             cout << "\t-h/-H/-?\tGive this help screen." << endl;
    813             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    814             cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
    815             cout << "\t-N\tGet non-convex-envelope." << endl;
    816             cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
    817             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
    818             cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
    819             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
    820             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
    821             cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
    822             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
    823             cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
    824             cout << "\t-v/-V\t\tGives version information." << endl;
    825             cout << "Note: config files must not begin with '-' !" << endl;
    826             delete(mol);
    827             delete(periode);
    828             return (1);
    829             break;
    830           case 'v':
    831           case 'V':
    832             cout << argv[0] << " " << VERSIONSTRING << endl;
    833             cout << "Build your own molecule position set." << endl;
    834             delete(mol);
    835             delete(periode);
    836             return (1);
    837             break;
    838           case 'e':
    839             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    840               cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
    841             } else {
    842               cout << "Using " << argv[argptr] << " as elements database." << endl;
    843               PathToDatabases = argv[argptr];
    844               argptr+=1;
    845             }
    846             break;
    847           case 'n':
    848             cout << "I won't parse trajectories." << endl;
    849             configuration.FastParsing = true;
    850             break;
    851           default:  // no match? Step on
    852             argptr++;
    853             break;
    854         }
    855       } else
    856         argptr++;
    857     } while (argptr < argc);
    858 
    859     // 2. Parse the element database
    860     if (periode->LoadPeriodentafel(PathToDatabases)) {
    861       cout << Verbose(0) << "Element list loaded successfully." << endl;
    862       //periode->Output((ofstream *)&cout);
    863     } else {
    864       cout << Verbose(0) << "Element list loading failed." << endl;
    865       return 1;
    866     }
    867 
    868     // 3. Find config file name and parse if possible
    869     if (argv[1][0] != '-') {
    870       cout << Verbose(0) << "Config file given." << endl;
    871       test.open(argv[1], ios::in);
    872       if (test == NULL) {
    873         //return (1);
    874         output.open(argv[1], ios::out);
    875         if (output == NULL) {
    876           cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
    877           config_present = absent;
    878         } else {
    879           cout << "Empty configuration file." << endl;
    880           ConfigFileName = argv[1];
    881           config_present = empty;
    882           output.close();
    883         }
    884       } else {
    885         test.close();
    886         ConfigFileName = argv[1];
    887         cout << Verbose(1) << "Specified config file found, parsing ... ";
    888         switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {
    889           case 1:
    890             cout << "new syntax." << endl;
    891             configuration.Load(ConfigFileName, periode, mol);
    892             config_present = present;
    893             break;
    894           case 0:
    895             cout << "old syntax." << endl;
    896             configuration.LoadOld(ConfigFileName, periode, mol);
    897             config_present = present;
    898             break;
    899           default:
    900             cout << "Unknown syntax or empty, yet present file." << endl;
    901             config_present = empty;
    902        }
    903       }
    904     } else
    905       config_present = absent;
    906 
    907     // 4. parse again through options, now for those depending on elements db and config presence
    908     argptr = 1;
    909     do {
    910       cout << "Current Command line argument: " << argv[argptr] << "." << endl;
    911       if (argv[argptr][0] == '-') {
    912         argptr++;
    913         if ((config_present == present) || (config_present == empty)) {
    914           switch(argv[argptr-1][1]) {
    915             case 'p':
    916               ExitFlag = 1;
    917               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    918                 ExitFlag = 255;
    919                 cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;
    920               } else {
    921                 SaveFlag = true;
    922                 cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
    923                 if (!mol->AddXYZFile(argv[argptr]))
    924                   cout << Verbose(2) << "File not found." << endl;
    925                 else {
    926                   cout << Verbose(2) << "File found and parsed." << endl;
    927                   config_present = present;
    928                 }
    929               }
    930               break;
    931             case 'a':
    932               ExitFlag = 1;
    933               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
    934                 ExitFlag = 255;
    935                 cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
    936               } else {
    937                 SaveFlag = true;
    938                 cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
    939                 first = new atom;
    940                 first->type = periode->FindElement(atoi(argv[argptr]));
    941                 if (first->type != NULL)
    942                   cout << Verbose(2) << "found element " << first->type->name << endl;
    943                 for (int i=NDIM;i--;)
    944                   first->x.x[i] = atof(argv[argptr+1+i]);
    945                 if (first->type != NULL) {
    946                   mol->AddAtom(first);  // add to molecule
    947                   if ((config_present == empty) && (mol->AtomCount != 0))
    948                     config_present = present;
    949                 } else
    950                   cerr << Verbose(1) << "Could not find the specified element." << endl;
    951                 argptr+=4;
    952               }
    953               break;
    954             default:   // no match? Don't step on (this is done in next switch's default)
    955               break;
    956           }
    957         }
    958         if (config_present == present) {
    959           switch(argv[argptr-1][1]) {
    960             case 'D':
    961               ExitFlag = 1;
    962               {
    963                 cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
    964                 MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
    965                 int *MinimumRingSize = new int[mol->AtomCount];
    966                 atom ***ListOfLocalAtoms = NULL;
    967                 int FragmentCounter = 0;
    968                 class StackClass<bond *> *BackEdgeStack = NULL;
    969                 class StackClass<bond *> *LocalBackEdgeStack = NULL;
    970                 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
    971                 mol->CreateListOfBondsPerAtom((ofstream *)&cout);
    972                 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
    973                 if (Subgraphs != NULL) {
    974                   Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
    975                   while (Subgraphs->next != NULL) {
    976                     Subgraphs = Subgraphs->next;
    977                     LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
    978                     Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
    979                     Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
    980                     delete(LocalBackEdgeStack);
    981                     delete(Subgraphs->previous);
    982                   }
    983                   delete(Subgraphs);
    984                   for (int i=0;i<FragmentCounter;i++)
    985                     Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");
    986                   Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");
    987                 }
    988                 delete(BackEdgeStack);
    989                 delete[](MinimumRingSize);
    990               }
    991               //argptr+=1;
    992               break;
    993             case 'E':
    994               ExitFlag = 1;
    995               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
    996                 ExitFlag = 255;
    997                 cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
    998               } else {
    999                 SaveFlag = true;
    1000                 cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
    1001                 first = mol->FindAtom(atoi(argv[argptr]));
    1002                 first->type = periode->FindElement(atoi(argv[argptr+1]));
    1003                 argptr+=2;
    1004               }
    1005               break;
    1006             case 'A':
    1007                 ExitFlag = 1;
    1008                 if ((argptr >= argc) || (argv[argptr][0] == '-')){
    1009                         ExitFlag =255;
    1010                         cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
    1011                 }
    1012                 else{
    1013                         cout << "Parsing bonds from " << argv[argptr] << "." << endl;
    1014                     ifstream *input = new ifstream(argv[argptr]);
    1015                         mol->CreateAdjacencyList2((ofstream *)&cout, input);
    1016                         input->close();
    1017                 }
    1018                 break;
    1019             case 'N':
    1020                 ExitFlag = 1;
    1021                 if ((argptr >= argc) || (argv[argptr][0] == '-')){
    1022                         ExitFlag = 255;
    1023                         cerr << "Not enough or invalid arguments given for non-convex envelope: -o <tecplot output file>" << endl;
    1024                         }
    1025                 else {
    1026                         cout << Verbose(0) << "Evaluating npn-convex envelope.";
    1027                         cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
    1028                         Find_non_convex_border((ofstream *)&cout, argv[argptr], mol);
    1029                         argptr+=1;
    1030                         }
    1031                 break;
    1032             case 'T':
    1033               ExitFlag = 1;
    1034               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1035                 ExitFlag = 255;
    1036                 cerr << "Not enough or invalid arguments given for storing temperature: -T <temperature file>" << endl;
    1037               } else {
    1038                 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
    1039                 ofstream *output = new ofstream(argv[argptr], ios::trunc);
    1040                 if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
    1041                   cout << Verbose(2) << "File could not be written." << endl;
    1042                 else
    1043                   cout << Verbose(2) << "File stored." << endl;
    1044                 output->close();
    1045                 delete(output);
    1046                 argptr+=1;
    1047               }
    1048               break;
    1049             case 'P':
    1050               ExitFlag = 1;
    1051               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1052                 ExitFlag = 255;
    1053                 cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
    1054               } else {
    1055                 SaveFlag = true;
    1056                 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
    1057                 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
    1058                   cout << Verbose(2) << "File not found." << endl;
    1059                 else
    1060                   cout << Verbose(2) << "File found and parsed." << endl;
    1061                 argptr+=1;
    1062               }
    1063               break;
    1064             case 't':
    1065               ExitFlag = 1;
    1066               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    1067                 ExitFlag = 255;
    1068                 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
    1069               } else {
    1070                 ExitFlag = 1;
    1071                 SaveFlag = true;
    1072                 cout << Verbose(1) << "Translating all ions to new origin." << endl;
    1073                 for (int i=NDIM;i--;)
    1074                   x.x[i] = atof(argv[argptr+i]);
    1075                 mol->Translate((const Vector *)&x);
    1076                 argptr+=3;
    1077               }
    1078               break;
    1079             case 's':
    1080               ExitFlag = 1;
    1081               if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
    1082                 ExitFlag = 255;
    1083                 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
    1084               } else {
    1085                 SaveFlag = true;
    1086                 j = -1;
    1087                 cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
    1088                 factor = new double[NDIM];
    1089                 factor[0] = atof(argv[argptr]);
    1090                 if ((argptr < argc) && (IsValidNumber(argv[argptr])))
    1091                   argptr++;
    1092                 factor[1] = atof(argv[argptr]);
    1093                 if ((argptr < argc) && (IsValidNumber(argv[argptr])))
    1094                   argptr++;
    1095                 factor[2] = atof(argv[argptr]);
    1096                 mol->Scale(&factor);
    1097                 for (int i=0;i<NDIM;i++) {
    1098                   j += i+1;
    1099                   x.x[i] = atof(argv[NDIM+i]);
    1100                   mol->cell_size[j]*=factor[i];
    1101                 }
    1102                 delete[](factor);
    1103                 argptr+=1;
    1104               }
    1105               break;
    1106             case 'b':
    1107               ExitFlag = 1;
    1108               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    1109                 ExitFlag = 255;
    1110                 cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl;
    1111               } else {
    1112                 SaveFlag = true;
    1113                 j = -1;
    1114                 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
    1115                 j=-1;
    1116                 for (int i=0;i<NDIM;i++) {
    1117                   j += i+1;
    1118                   x.x[i] = atof(argv[argptr++]);
    1119                   mol->cell_size[j] += x.x[i]*2.;
    1120                 }
    1121                 // center
    1122                 mol->CenterInBox((ofstream *)&cout, &x);
    1123                 // update Box of atoms by boundary
    1124                 mol->SetBoxDimension(&x);
    1125               }
    1126               break;
    1127             case 'c':
    1128               ExitFlag = 1;
    1129               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    1130                 ExitFlag = 255;
    1131                 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
    1132               } else {
    1133                 SaveFlag = true;
    1134                 j = -1;
    1135                 cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
    1136                 // make every coordinate positive
    1137                 mol->CenterEdge((ofstream *)&cout, &x);
    1138                 // update Box of atoms by boundary
    1139                 mol->SetBoxDimension(&x);
    1140                 // translate each coordinate by boundary
    1141                 j=-1;
    1142                 for (int i=0;i<NDIM;i++) {
    1143                   j += i+1;
    1144                   x.x[i] = atof(argv[argptr++]);
    1145                   mol->cell_size[j] += x.x[i]*2.;
    1146                 }
    1147                 mol->Translate((const Vector *)&x);
    1148               }
    1149               break;
    1150             case 'O':
    1151               ExitFlag = 1;
    1152               SaveFlag = true;
    1153               cout << Verbose(1) << "Centering atoms in origin." << endl;
    1154               mol->CenterOrigin((ofstream *)&cout, &x);
    1155               mol->SetBoxDimension(&x);
    1156               break;
    1157             case 'r':
    1158               ExitFlag = 1;
    1159               SaveFlag = true;
    1160               cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
    1161               break;
    1162             case 'F':
    1163             case 'f':
    1164               ExitFlag = 1;
    1165               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
    1166                 ExitFlag = 255;
    1167                 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
    1168               } else {
    1169                 cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
    1170                 cout << Verbose(0) << "Creating connection matrix..." << endl;
    1171                 start = clock();
    1172                 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
    1173                 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
    1174                 if (mol->first->next != mol->last) {
    1175                   ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
    1176                 }
    1177                 end = clock();
    1178                 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    1179                 argptr+=2;
    1180               }
    1181               break;
    1182             case 'm':
    1183               ExitFlag = 1;
    1184               j = atoi(argv[argptr++]);
    1185               if ((j<0) || (j>1)) {
    1186                 cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
    1187                 j = 0;
    1188               }
    1189               if (j) {
    1190                 SaveFlag = true;
    1191                 cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
    1192               } else
    1193                 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
    1194               mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
    1195               break;
    1196             case 'o':
    1197               ExitFlag = 1;
    1198               if ((argptr >= argc) || (argv[argptr][0] == '-')){
    1199                 ExitFlag = 255;
    1200                 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
    1201               } else {
    1202                 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
    1203                 ofstream *output = new ofstream(argv[argptr], ios::trunc);
    1204                 //$$$
    1205                 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
    1206                 VolumeOfConvexEnvelope((ofstream *)&cout, output, &configuration, NULL, mol);
    1207                 output->close();
    1208                 delete(output);
    1209                 argptr+=1;
    1210               }
    1211               break;
    1212             case 'U':
    1213               ExitFlag = 1;
    1214               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
    1215                 ExitFlag = 255;
    1216                 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
    1217                 volume = -1; // for case 'u': don't print error again
    1218               } else {
    1219                 volume = atof(argv[argptr++]);
    1220                 cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
    1221               }
    1222             case 'u':
    1223               ExitFlag = 1;
    1224               if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
    1225                 if (volume != -1)
    1226                   ExitFlag = 255;
    1227                   cerr << "Not enough arguments given for suspension: -u <density>" << endl;
    1228               } else {
    1229                 double density;
    1230                 SaveFlag = true;
    1231                 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
    1232                 density = atof(argv[argptr++]);
    1233                 if (density < 1.0) {
    1234                   cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
    1235                   density = 1.3;
    1236                 }
    1237 //                for(int i=0;i<NDIM;i++) {
    1238 //                  repetition[i] = atoi(argv[argptr++]);
    1239 //                  if (repetition[i] < 1)
    1240 //                    cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
    1241 //                  repetition[i] = 1;
    1242 //                }
    1243                 PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);  // if volume == 0, will calculate from ConvexEnvelope
    1244               }
    1245               break;
    1246             case 'd':
    1247               ExitFlag = 1;
    1248               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    1249                 ExitFlag = 255;
    1250                 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
    1251               } else {
    1252                 SaveFlag = true;
    1253                 for (int axis = 1; axis <= NDIM; axis++) {
    1254                   int faktor = atoi(argv[argptr++]);
    1255                   int count;
    1256                   element ** Elements;
    1257                   Vector ** vectors;
    1258                   if (faktor < 1) {
    1259                     cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
    1260                     faktor = 1;
    1261                   }
    1262                   mol->CountAtoms((ofstream *)&cout);  // recount atoms
    1263                   if (mol->AtomCount != 0) {  // if there is more than none
    1264                     count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
    1265                     Elements = new element *[count];
    1266                     vectors = new Vector *[count];
    1267                     j = 0;
    1268                     first = mol->start;
    1269                     while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
    1270                       first = first->next;
    1271                       Elements[j] = first->type;
    1272                       vectors[j] = &first->x;
    1273                       j++;
    1274                     }
    1275                     if (count != j)
    1276                       cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
    1277                     x.Zero();
    1278                     y.Zero();
    1279                     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
    1280                     for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    1281                       x.AddVector(&y); // per factor one cell width further
    1282                       for (int k=count;k--;) { // go through every atom of the original cell
    1283                         first = new atom(); // create a new body
    1284                         first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    1285                         first->x.AddVector(&x);      // translate the coordinates
    1286                         first->type = Elements[k];  // insert original element
    1287                         mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
    1288                       }
    1289                     }
    1290                     // free memory
    1291                     delete[](Elements);
    1292                     delete[](vectors);
    1293                     // correct cell size
    1294                     if (axis < 0) { // if sign was negative, we have to translate everything
    1295                       x.Zero();
    1296                       x.AddVector(&y);
    1297                       x.Scale(-(faktor-1));
    1298                       mol->Translate(&x);
    1299                     }
    1300                     mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    1301                   }
    1302                 }
    1303               }
    1304               break;
    1305             default:  // no match? Step on
    1306               if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
    1307                 argptr++;
    1308               break;
    1309           }
    1310         }
    1311       } else argptr++;
    1312     } while (argptr < argc);
    1313     if (SaveFlag)
    1314       SaveConfig(ConfigFileName, &configuration, periode, mol);
    1315     if ((ExitFlag >= 1)) {
    1316       delete(mol);
    1317       delete(periode);
    1318       return (ExitFlag);
    1319     }
    1320   } else {  // no arguments, hence scan the elements db
    1321     if (periode->LoadPeriodentafel(PathToDatabases))
    1322       cout << Verbose(0) << "Element list loaded successfully." << endl;
    1323     else
    1324       cout << Verbose(0) << "Element list loading failed." << endl;
    1325     configuration.RetrieveConfigPathAndName("main_pcp_linux");
    1326   }
    1327   return(0);
     777        Vector x,y,z,n; // coordinates for absolute point in cell volume
     778        double *factor; // unit factor if desired
     779        ifstream test;
     780        ofstream output;
     781        string line;
     782        atom *first;
     783        bool SaveFlag = false;
     784        int ExitFlag = 0;
     785        int j;
     786        double volume = 0.;
     787        enum ConfigStatus config_present = absent;
     788        clock_t start,end;
     789        int argptr;
     790        PathToDatabases = LocalPath;
     791
     792        if (argc > 1) { // config file specified as option
     793                // 1. : Parse options that just set variables or print help
     794                argptr = 1;
     795                do {
     796                        if (argv[argptr][0] == '-') {
     797                                cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
     798                                argptr++;
     799                                switch(argv[argptr-1][1]) {
     800                                        case 'h':
     801                                        case 'H':
     802                                        case '?':
     803                                                cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
     804                                                cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
     805                                                cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
     806                                                cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
     807                                                cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
     808                                                cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
     809                                                cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
     810                                                cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
     811                                                cout << "\t-O\tCenter atoms in origin." << endl;
     812                                                cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
     813                                                cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
     814                                                cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
     815                                                cout << "\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;
     816                                                cout << "\t-h/-H/-?\tGive this help screen." << endl;
     817                                                cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
     818                                                cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
     819                                                cout << "\t-N\tGet non-convex-envelope." << endl;
     820                                                cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
     821                                                cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
     822                                                cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
     823                                                cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
     824                                                cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
     825                                                cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
     826                                                cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
     827                                                cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
     828                                                cout << "\t-v/-V\t\tGives version information." << endl;
     829                                                cout << "Note: config files must not begin with '-' !" << endl;
     830                                                delete(mol);
     831                                                delete(periode);
     832                                                return (1);
     833                                                break;
     834                                        case 'v':
     835                                        case 'V':
     836                                                cout << argv[0] << " " << VERSIONSTRING << endl;
     837                                                cout << "Build your own molecule position set." << endl;
     838                                                delete(mol);
     839                                                delete(periode);
     840                                                return (1);
     841                                                break;
     842                                        case 'e':
     843                                                if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     844                                                        cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
     845                                                } else {
     846                                                        cout << "Using " << argv[argptr] << " as elements database." << endl;
     847                                                        PathToDatabases = argv[argptr];
     848                                                        argptr+=1;
     849                                                }
     850                                                break;
     851                                        case 'n':
     852                                                cout << "I won't parse trajectories." << endl;
     853                                                configuration.FastParsing = true;
     854                                                break;
     855                                        default:        // no match? Step on
     856                                                argptr++;
     857                                                break;
     858                                }
     859                        } else
     860                                argptr++;
     861                } while (argptr < argc);
     862
     863                // 2. Parse the element database
     864                if (periode->LoadPeriodentafel(PathToDatabases)) {
     865                        cout << Verbose(0) << "Element list loaded successfully." << endl;
     866                        //periode->Output((ofstream *)&cout);
     867                } else {
     868                        cout << Verbose(0) << "Element list loading failed." << endl;
     869                        return 1;
     870                }
     871                // 3. Find config file name and parse if possible
     872                if (argv[1][0] != '-') {
     873                        cout << Verbose(0) << "Config file given." << endl;
     874                        test.open(argv[1], ios::in);
     875                        if (test == NULL) {
     876                                //return (1);
     877                                output.open(argv[1], ios::out);
     878                                if (output == NULL) {
     879                                        cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
     880                                        config_present = absent;
     881                                } else {
     882                                        cout << "Empty configuration file." << endl;
     883                                        ConfigFileName = argv[1];
     884                                        config_present = empty;
     885                                        output.close();
     886                                }
     887                        } else {
     888                                test.close();
     889                                ConfigFileName = argv[1];
     890                                cout << Verbose(1) << "Specified config file found, parsing ... ";
     891                                switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {
     892                                        case 1:
     893                                                cout << "new syntax." << endl;
     894                                                configuration.Load(ConfigFileName, periode, mol);
     895                                                config_present = present;
     896                                                break;
     897                                        case 0:
     898                                                cout << "old syntax." << endl;
     899                                                configuration.LoadOld(ConfigFileName, periode, mol);
     900                                                config_present = present;
     901                                                break;
     902                                        default:
     903                                                cout << "Unknown syntax or empty, yet present file." << endl;
     904                                                config_present = empty;
     905                         }
     906                        }
     907                } else
     908                        config_present = absent;
     909                // 4. parse again through options, now for those depending on elements db and config presence
     910                argptr = 1;
     911                do {
     912                        cout << "Current Command line argument: " << argv[argptr] << "." << endl;
     913                        if (argv[argptr][0] == '-') {
     914                                argptr++;
     915                                if ((config_present == present) || (config_present == empty)) {
     916                                        switch(argv[argptr-1][1]) {
     917                                                case 'p':
     918                                                        ExitFlag = 1;
     919                                                        if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     920                                                                ExitFlag = 255;
     921                                                                cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;
     922                                                        } else {
     923                                                                SaveFlag = true;
     924                                                                cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
     925                                                                if (!mol->AddXYZFile(argv[argptr]))
     926                                                                        cout << Verbose(2) << "File not found." << endl;
     927                                                                else {
     928                                                                        cout << Verbose(2) << "File found and parsed." << endl;
     929                                                                        config_present = present;
     930                                                                }
     931                                                        }
     932                                                        break;
     933                                                case 'a':
     934                                                        ExitFlag = 1;
     935                                                        if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
     936                                                                ExitFlag = 255;
     937                                                                cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
     938                                                        } else {
     939                                                                SaveFlag = true;
     940                                                                cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
     941                                                                first = new atom;
     942                                                                first->type = periode->FindElement(atoi(argv[argptr]));
     943                                                                if (first->type != NULL)
     944                                                                        cout << Verbose(2) << "found element " << first->type->name << endl;
     945                                                                for (int i=NDIM;i--;)
     946                                                                        first->x.x[i] = atof(argv[argptr+1+i]);
     947                                                                if (first->type != NULL) {
     948                                                                        mol->AddAtom(first);    // add to molecule
     949                                                                        if ((config_present == empty) && (mol->AtomCount != 0))
     950                                                                                config_present = present;
     951                                                                } else
     952                                                                        cerr << Verbose(1) << "Could not find the specified element." << endl;
     953                                                                argptr+=4;
     954                                                        }
     955                                                        break;
     956                                                default:         // no match? Don't step on (this is done in next switch's default)
     957                                                        break;
     958                                        }
     959                                }
     960                                if (config_present == present) {
     961                                        switch(argv[argptr-1][1]) {
     962                                                case 'D':
     963                                                        ExitFlag = 1;
     964                                                        {
     965                                                                cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
     966                                                                MoleculeLeafClass *Subgraphs = NULL;                    // list of subgraphs from DFS analysis
     967                                                                int *MinimumRingSize = new int[mol->AtomCount];
     968                                                                atom ***ListOfLocalAtoms = NULL;
     969                                                                int FragmentCounter = 0;
     970                                                                class StackClass<bond *> *BackEdgeStack = NULL;
     971                                                                class StackClass<bond *> *LocalBackEdgeStack = NULL;
     972                                                                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
     973                                                                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
     974                                                                Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
     975                                                                if (Subgraphs != NULL) {
     976                                                                        Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false);        // we want to keep the created ListOfLocalAtoms
     977                                                                        while (Subgraphs->next != NULL) {
     978                                                                                Subgraphs = Subgraphs->next;
     979                                                                                LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
     980                                                                                Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
     981                                                                                Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
     982                                                                                delete(LocalBackEdgeStack);
     983                                                                                delete(Subgraphs->previous);
     984                                                                        }
     985                                                                        delete(Subgraphs);
     986                                                                        for (int i=0;i<FragmentCounter;i++)
     987                                                                                Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");
     988                                                                        Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");
     989                                                                }
     990                                                                delete(BackEdgeStack);
     991                                                                delete[](MinimumRingSize);
     992                                                        }
     993                                                        //argptr+=1;
     994                                                        break;
     995                                                case 'E':
     996                                                        ExitFlag = 1;
     997                                                        if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
     998                                                                ExitFlag = 255;
     999                                                                cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
     1000                                                        } else {
     1001                                                                SaveFlag = true;
     1002                                                                cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
     1003                                                                first = mol->FindAtom(atoi(argv[argptr]));
     1004                                                                first->type = periode->FindElement(atoi(argv[argptr+1]));
     1005                                                                argptr+=2;
     1006                                                        }
     1007                                                        break;
     1008                                                case 'A':
     1009                                                        ExitFlag = 1;
     1010                                                        if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1011                                                                ExitFlag =255;
     1012                                                                cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
     1013                                                        } else {
     1014                                                                cout << "Parsing bonds from " << argv[argptr] << "." << endl;
     1015                                                                ifstream *input = new ifstream(argv[argptr]);
     1016                                                                mol->CreateAdjacencyList2((ofstream *)&cout, input);
     1017                                                                input->close();
     1018                                                                argptr+=1;
     1019                                                        }
     1020                                                        break;
     1021                                                case 'N':
     1022                                                        ExitFlag = 1;
     1023                                                        if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
     1024                                                                ExitFlag = 255;
     1025                                                                cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
     1026                                                        } else {
     1027                                                                class Tesselation T;
     1028                                                                int N = 15;
     1029                                                                int number = 100;
     1030                                                                string filename(argv[argptr+1]);
     1031                                                                filename.append(".csv");
     1032                                                                cout << Verbose(0) << "Evaluating non-convex envelope.";
     1033                                                                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
     1034                                                                LinkedCell LCList(mol, atof(argv[argptr]));             // \NOTE not twice the radius??
     1035                                                                Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr]));
     1036                                                                FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
     1037                                                                argptr+=2;
     1038                                                        }
     1039                                                        break;
     1040                                                case 'T':
     1041                                                        ExitFlag = 1;
     1042                                                        if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1043                                                                ExitFlag = 255;
     1044                                                                cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
     1045                                                        } else {
     1046                                                                cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
     1047                                                                ofstream *output = new ofstream(argv[argptr], ios::trunc);
     1048                                                                if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
     1049                                                                        cout << Verbose(2) << "File could not be written." << endl;
     1050                                                                else
     1051                                                                        cout << Verbose(2) << "File stored." << endl;
     1052                                                                output->close();
     1053                                                                delete(output);
     1054                                                                argptr+=1;
     1055                                                        }
     1056                                                        break;
     1057                                                case 'P':
     1058                                                        ExitFlag = 1;
     1059                                                        if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1060                                                                ExitFlag = 255;
     1061                                                                cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
     1062                                                        } else {
     1063                                                                SaveFlag = true;
     1064                                                                cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
     1065                                                                if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
     1066                                                                        cout << Verbose(2) << "File not found." << endl;
     1067                                                                else
     1068                                                                        cout << Verbose(2) << "File found and parsed." << endl;
     1069                                                                argptr+=1;
     1070                                                        }
     1071                                                        break;
     1072                                                case 't':
     1073                                                        ExitFlag = 1;
     1074                                                        if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
     1075                                                                ExitFlag = 255;
     1076                                                                cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
     1077                                                        } else {
     1078                                                                ExitFlag = 1;
     1079                                                                SaveFlag = true;
     1080                                                                cout << Verbose(1) << "Translating all ions to new origin." << endl;
     1081                                                                for (int i=NDIM;i--;)
     1082                                                                        x.x[i] = atof(argv[argptr+i]);
     1083                                                                mol->Translate((const Vector *)&x);
     1084                                                                argptr+=3;
     1085                                                        }
     1086                                                        break;
     1087                                                case 's':
     1088                                                        ExitFlag = 1;
     1089                                                        if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
     1090                                                                ExitFlag = 255;
     1091                                                                cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
     1092                                                        } else {
     1093                                                                SaveFlag = true;
     1094                                                                j = -1;
     1095                                                                cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
     1096                                                                factor = new double[NDIM];
     1097                                                                factor[0] = atof(argv[argptr]);
     1098                                                                if ((argptr < argc) && (IsValidNumber(argv[argptr])))
     1099                                                                        argptr++;
     1100                                                                factor[1] = atof(argv[argptr]);
     1101                                                                if ((argptr < argc) && (IsValidNumber(argv[argptr])))
     1102                                                                        argptr++;
     1103                                                                factor[2] = atof(argv[argptr]);
     1104                                                                mol->Scale(&factor);
     1105                                                                for (int i=0;i<NDIM;i++) {
     1106                                                                        j += i+1;
     1107                                                                        x.x[i] = atof(argv[NDIM+i]);
     1108                                                                        mol->cell_size[j]*=factor[i];
     1109                                                                }
     1110                                                                delete[](factor);
     1111                                                                argptr+=1;
     1112                                                        }
     1113                                                        break;
     1114                                                case 'b':
     1115                                                        ExitFlag = 1;
     1116                                                        if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
     1117                                                                ExitFlag = 255;
     1118                                                                cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl;
     1119                                                        } else {
     1120                                                                SaveFlag = true;
     1121                                                                j = -1;
     1122                                                                cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     1123                                                                j=-1;
     1124                                                                for (int i=0;i<NDIM;i++) {
     1125                                                                        j += i+1;
     1126                                                                        x.x[i] = atof(argv[argptr++]);
     1127                                                                        mol->cell_size[j] += x.x[i]*2.;
     1128                                                                }
     1129                                                                // center
     1130                                                                mol->CenterInBox((ofstream *)&cout, &x);
     1131                                                                // update Box of atoms by boundary
     1132                                                                mol->SetBoxDimension(&x);
     1133                                                        }
     1134                                                        break;
     1135                                                case 'c':
     1136                                                        ExitFlag = 1;
     1137                                                        if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
     1138                                                                ExitFlag = 255;
     1139                                                                cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
     1140                                                        } else {
     1141                                                                SaveFlag = true;
     1142                                                                j = -1;
     1143                                                                cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
     1144                                                                // make every coordinate positive
     1145                                                                mol->CenterEdge((ofstream *)&cout, &x);
     1146                                                                // update Box of atoms by boundary
     1147                                                                mol->SetBoxDimension(&x);
     1148                                                                // translate each coordinate by boundary
     1149                                                                j=-1;
     1150                                                                for (int i=0;i<NDIM;i++) {
     1151                                                                        j += i+1;
     1152                                                                        x.x[i] = atof(argv[argptr++]);
     1153                                                                        mol->cell_size[j] += x.x[i]*2.;
     1154                                                                }
     1155                                                                mol->Translate((const Vector *)&x);
     1156                                                        }
     1157                                                        break;
     1158                                                case 'O':
     1159                                                        ExitFlag = 1;
     1160                                                        SaveFlag = true;
     1161                                                        cout << Verbose(1) << "Centering atoms in origin." << endl;
     1162                                                        mol->CenterOrigin((ofstream *)&cout, &x);
     1163                                                        mol->SetBoxDimension(&x);
     1164                                                        break;
     1165                                                case 'r':
     1166                                                        ExitFlag = 1;
     1167                                                        SaveFlag = true;
     1168                                                        cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
     1169                                                        break;
     1170                                                case 'F':
     1171                                                case 'f':
     1172                                                        ExitFlag = 1;
     1173                                                        if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
     1174                                                                ExitFlag = 255;
     1175                                                                cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
     1176                                                        } else {
     1177                                                                cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
     1178                                                                cout << Verbose(0) << "Creating connection matrix..." << endl;
     1179                                                                start = clock();
     1180                                                                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
     1181                                                                cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
     1182                                                                if (mol->first->next != mol->last) {
     1183                                                                        ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
     1184                                                                }
     1185                                                                end = clock();
     1186                                                                cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     1187                                                                argptr+=2;
     1188                                                        }
     1189                                                        break;
     1190                                                case 'm':
     1191                                                        ExitFlag = 1;
     1192                                                        j = atoi(argv[argptr++]);
     1193                                                        if ((j<0) || (j>1)) {
     1194                                                                cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
     1195                                                                j = 0;
     1196                                                        }
     1197                                                        if (j) {
     1198                                                                SaveFlag = true;
     1199                                                                cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
     1200                                                        } else
     1201                                                                cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
     1202                                                        mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
     1203                                                        break;
     1204                                                case 'o':
     1205                                                        ExitFlag = 1;
     1206                                                        if ((argptr >= argc) || (argv[argptr][0] == '-')){
     1207                                                                ExitFlag = 255;
     1208                                                                cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
     1209                                                        } else {
     1210                                                                cout << Verbose(0) << "Evaluating volume of the convex envelope.";
     1211                                                                cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
     1212                                                                VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol);
     1213                                                                argptr+=1;
     1214                                                        }
     1215                                                        break;
     1216                                                case 'U':
     1217                                                        ExitFlag = 1;
     1218                                                        if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
     1219                                                                ExitFlag = 255;
     1220                                                                cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
     1221                                                                volume = -1; // for case 'u': don't print error again
     1222                                                        } else {
     1223                                                                volume = atof(argv[argptr++]);
     1224                                                                cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
     1225                                                        }
     1226                                                case 'u':
     1227                                                        ExitFlag = 1;
     1228                                                        if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
     1229                                                                if (volume != -1)
     1230                                                                        ExitFlag = 255;
     1231                                                                        cerr << "Not enough arguments given for suspension: -u <density>" << endl;
     1232                                                        } else {
     1233                                                                double density;
     1234                                                                SaveFlag = true;
     1235                                                                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
     1236                                                                density = atof(argv[argptr++]);
     1237                                                                if (density < 1.0) {
     1238                                                                        cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
     1239                                                                        density = 1.3;
     1240                                                                }
     1241//                                                              for(int i=0;i<NDIM;i++) {
     1242//                                                                      repetition[i] = atoi(argv[argptr++]);
     1243//                                                                      if (repetition[i] < 1)
     1244//                                                                              cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
     1245//                                                                      repetition[i] = 1;
     1246//                                                              }
     1247                                                                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);        // if volume == 0, will calculate from ConvexEnvelope
     1248                                                        }
     1249                                                        break;
     1250                                                case 'd':
     1251                                                        ExitFlag = 1;
     1252                                                        if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
     1253                                                                ExitFlag = 255;
     1254                                                                cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
     1255                                                        } else {
     1256                                                                SaveFlag = true;
     1257                                                                for (int axis = 1; axis <= NDIM; axis++) {
     1258                                                                        int faktor = atoi(argv[argptr++]);
     1259                                                                        int count;
     1260                                                                        element ** Elements;
     1261                                                                        Vector ** vectors;
     1262                                                                        if (faktor < 1) {
     1263                                                                                cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
     1264                                                                                faktor = 1;
     1265                                                                        }
     1266                                                                        mol->CountAtoms((ofstream *)&cout);     // recount atoms
     1267                                                                        if (mol->AtomCount != 0) {      // if there is more than none
     1268                                                                                count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
     1269                                                                                Elements = new element *[count];
     1270                                                                                vectors = new Vector *[count];
     1271                                                                                j = 0;
     1272                                                                                first = mol->start;
     1273                                                                                while (first->next != mol->end) {       // make a list of all atoms with coordinates and element
     1274                                                                                        first = first->next;
     1275                                                                                        Elements[j] = first->type;
     1276                                                                                        vectors[j] = &first->x;
     1277                                                                                        j++;
     1278                                                                                }
     1279                                                                                if (count != j)
     1280                                                                                        cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
     1281                                                                                x.Zero();
     1282                                                                                y.Zero();
     1283                                                                                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
     1284                                                                                for (int i=1;i<faktor;i++) {    // then add this list with respective translation factor times
     1285                                                                                        x.AddVector(&y); // per factor one cell width further
     1286                                                                                        for (int k=count;k--;) { // go through every atom of the original cell
     1287                                                                                                first = new atom(); // create a new body
     1288                                                                                                first->x.CopyVector(vectors[k]);        // use coordinate of original atom
     1289                                                                                                first->x.AddVector(&x);                 // translate the coordinates
     1290                                                                                                first->type = Elements[k];      // insert original element
     1291                                                                                                mol->AddAtom(first);                            // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
     1292                                                                                        }
     1293                                                                                }
     1294                                                                                // free memory
     1295                                                                                delete[](Elements);
     1296                                                                                delete[](vectors);
     1297                                                                                // correct cell size
     1298                                                                                if (axis < 0) { // if sign was negative, we have to translate everything
     1299                                                                                        x.Zero();
     1300                                                                                        x.AddVector(&y);
     1301                                                                                        x.Scale(-(faktor-1));
     1302                                                                                        mol->Translate(&x);
     1303                                                                                }
     1304                                                                                mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     1305                                                                        }
     1306                                                                }
     1307                                                        }
     1308                                                        break;
     1309                                                default:        // no match? Step on
     1310                                                        if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
     1311                                                                argptr++;
     1312                                                        break;
     1313                                        }
     1314                                }
     1315                        } else argptr++;
     1316                } while (argptr < argc);
     1317                if (SaveFlag)
     1318                        SaveConfig(ConfigFileName, &configuration, periode, mol);
     1319                if ((ExitFlag >= 1)) {
     1320                        delete(mol);
     1321                        delete(periode);
     1322                        return (ExitFlag);
     1323                }
     1324        } else {        // no arguments, hence scan the elements db
     1325                if (periode->LoadPeriodentafel(PathToDatabases))
     1326                        cout << Verbose(0) << "Element list loaded successfully." << endl;
     1327                else
     1328                        cout << Verbose(0) << "Element list loading failed." << endl;
     1329                configuration.RetrieveConfigPathAndName("main_pcp_linux");
     1330        }
     1331        return(0);
    13281332};
    13291333
    … …  
    13321336int main(int argc, char **argv)
    13331337{
    1334   periodentafel *periode = new periodentafel; // and a period table of all elements
    1335   molecule *mol = new molecule(periode);    // first we need an empty molecule
    1336   config configuration;
    1337   double tmp1;
    1338   double bond, min_bond;
    1339   atom *first, *second;
    1340   char choice;  // menu choice char
    1341   Vector x,y,z,n;  // coordinates for absolute point in cell volume
     1338        periodentafel *periode = new periodentafel; // and a period table of all elements
     1339        molecule *mol = new molecule(periode);          // first we need an empty molecule
     1340        config configuration;
     1341        double tmp1;
     1342        double bond, min_bond;
     1343        atom *first, *second;
     1344        char choice;    // menu choice char
     1345        Vector x,y,z,n; // coordinates for absolute point in cell volume
    13421346        double *factor; // unit factor if desired
    1343   bool valid; // flag if input was valid or not
    1344   ifstream test;
    1345   ofstream output;
    1346   string line;
    1347   char filename[MAXSTRINGSIZE];
    1348   char *ConfigFileName = NULL;
    1349   char *ElementsFileName = NULL;
    1350   int Z;
    1351   int j, axis, count, faktor;
    1352   clock_t start,end;
    1353 //  int *MinimumRingSize = NULL;
    1354   MoleculeLeafClass *Subgraphs = NULL;
    1355 //  class StackClass<bond *> *BackEdgeStack = NULL;
    1356   element **Elements;
    1357   Vector **vectors;
    1358 
    1359   // =========================== PARSE COMMAND LINE OPTIONS ====================================
    1360   j = ParseCommandLineOptions(argc, argv, mol, periode, configuration, ConfigFileName, ElementsFileName);
    1361   if (j == 1) return 0; // just for -v and -h options
    1362   if (j) return j;  // something went wrong
    1363 
    1364   // General stuff
    1365   if (mol->cell_size[0] == 0.) {
    1366     cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
    1367     for (int i=0;i<6;i++) {
    1368       cout << Verbose(1) << "Cell size" << i << ": ";
    1369       cin >> mol->cell_size[i];
    1370     }
    1371   }
    1372 
    1373   // =========================== START INTERACTIVE SESSION ====================================
    1374 
    1375   // now the main construction loop
    1376   cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
    1377   do {
    1378     cout << Verbose(0) << endl << endl;
    1379     cout << Verbose(0) << "============Element list=======================" << endl;
    1380     mol->Checkout((ofstream *)&cout);
    1381     cout << Verbose(0) << "============Atom list==========================" << endl;
    1382     mol->Output((ofstream *)&cout);
    1383     cout << Verbose(0) << "============Menu===============================" << endl;
    1384     cout << Verbose(0) << "a - add an atom" << endl;
    1385     cout << Verbose(0) << "r - remove an atom" << endl;
    1386     cout << Verbose(0) << "b - scale a bond between atoms" << endl;
    1387     cout << Verbose(0) << "u - change an atoms element" << endl;
    1388     cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
    1389     cout << Verbose(0) << "-----------------------------------------------" << endl;
    1390     cout << Verbose(0) << "p - Parse xyz file" << endl;
    1391     cout << Verbose(0) << "e - edit the current configuration" << endl;
    1392     cout << Verbose(0) << "o - create connection matrix" << endl;
    1393     cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
    1394     cout << Verbose(0) << "-----------------------------------------------" << endl;
    1395     cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
    1396     cout << Verbose(0) << "i - realign molecule" << endl;
    1397     cout << Verbose(0) << "m - mirror all molecules" << endl;
    1398     cout << Verbose(0) << "t - translate molecule by vector" << endl;
    1399     cout << Verbose(0) << "c - scale by unit transformation" << endl;
    1400     cout << Verbose(0) << "g - center atoms in box" << endl;
    1401     cout << Verbose(0) << "-----------------------------------------------" << endl;
    1402     cout << Verbose(0) << "s - save current setup to config file" << endl;
    1403     cout << Verbose(0) << "T - call the current test routine" << endl;
    1404     cout << Verbose(0) << "q - quit" << endl;
    1405     cout << Verbose(0) << "===============================================" << endl;
    1406     cout << Verbose(0) << "Input: ";
    1407     cin >> choice;
    1408 
    1409     switch (choice) {
    1410       default:
    1411       case 'a': // add atom
    1412         AddAtoms(periode, mol);
    1413         choice = 'a';
    1414         break;
    1415 
    1416       case 'b': // scale a bond
    1417         cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
    1418         first = mol->AskAtom("Enter first (fixed) atom: ");
    1419         second = mol->AskAtom("Enter second (shifting) atom: ");
    1420         min_bond = 0.;
    1421         for (int i=NDIM;i--;)
    1422           min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
    1423         min_bond = sqrt(min_bond);
    1424         cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
    1425         cout << Verbose(0) << "Enter new bond length [a.u.]: ";
    1426         cin >> bond;
    1427         for (int i=NDIM;i--;) {
    1428           second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
    1429         }
    1430         //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
    1431         //second->Output(second->type->No, 1, (ofstream *)&cout);
    1432         break;
    1433 
    1434       case 'c': // unit scaling of the metric
    1435       cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
    1436       cout << Verbose(0) << "Enter three factors: ";
    1437       factor = new double[NDIM];
    1438       cin >> factor[0];
    1439       cin >> factor[1];
    1440       cin >> factor[2];
    1441       valid = true;
    1442       mol->Scale(&factor);
    1443       delete[](factor);
    1444       break;
    1445 
    1446       case 'd': // duplicate the periodic cell along a given axis, given times
    1447         cout << Verbose(0) << "State the axis [(+-)123]: ";
    1448         cin >> axis;
    1449         cout << Verbose(0) << "State the factor: ";
    1450         cin >> faktor;
    1451 
    1452         mol->CountAtoms((ofstream *)&cout);  // recount atoms
    1453         if (mol->AtomCount != 0) {  // if there is more than none
    1454           count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
    1455           Elements = new element *[count];
    1456           vectors = new Vector *[count];
    1457           j = 0;
    1458           first = mol->start;
    1459           while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
    1460             first = first->next;
    1461             Elements[j] = first->type;
    1462             vectors[j] = &first->x;
    1463             j++;
    1464           }
    1465           if (count != j)
    1466             cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
    1467           x.Zero();
    1468           y.Zero();
    1469           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
    1470           for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    1471             x.AddVector(&y); // per factor one cell width further
    1472             for (int k=count;k--;) { // go through every atom of the original cell
    1473               first = new atom(); // create a new body
    1474               first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    1475               first->x.AddVector(&x);      // translate the coordinates
    1476               first->type = Elements[k];  // insert original element
    1477               mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
    1478             }
    1479           }
    1480           if (mol->first->next != mol->last) // if connect matrix is present already, redo it
    1481             mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration.GetIsAngstroem());
    1482           // free memory
    1483           delete[](Elements);
    1484           delete[](vectors);
    1485           // correct cell size
    1486           if (axis < 0) { // if sign was negative, we have to translate everything
    1487             x.Zero();
    1488             x.AddVector(&y);
    1489             x.Scale(-(faktor-1));
    1490             mol->Translate(&x);
    1491           }
    1492           mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    1493         }
    1494         break;
    1495 
    1496       case 'e': // edit each field of the configuration
    1497       configuration.Edit(mol);
    1498       break;
    1499 
    1500       case 'f':
    1501         FragmentAtoms(mol, &configuration);
    1502         break;
    1503 
    1504       case 'g': // center the atoms
    1505         CenterAtoms(mol);
    1506         break;
    1507 
    1508       case 'i': // align all atoms
    1509         AlignAtoms(periode, mol);
    1510         break;
    1511 
    1512       case 'l': // measure distances or angles
    1513         MeasureAtoms(periode, mol, &configuration);
    1514         break;
    1515 
    1516       case 'm': // mirror atoms along a given axis
    1517         MirrorAtoms(mol);
    1518         break;
    1519 
    1520       case 'o': // create the connection matrix
    1521         {
    1522           cout << Verbose(0) << "What's the maximum bond distance: ";
    1523           cin >> tmp1;
    1524           start = clock();
    1525           mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
    1526           mol->CreateListOfBondsPerAtom((ofstream *)&cout);
    1527 //          Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
    1528 //          while (Subgraphs->next != NULL) {
    1529 //            Subgraphs = Subgraphs->next;
    1530 //            Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
    1531 //            delete(Subgraphs->previous);
    1532 //          }
    1533 //          delete(Subgraphs);    // we don't need the list here, so free everything
    1534 //          delete[](MinimumRingSize);
    1535 //          Subgraphs = NULL;
    1536           end = clock();
    1537           cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    1538         }
    1539         break;
    1540 
    1541       case 'p': // parse and XYZ file
    1542         cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
    1543         do {
    1544           cout << Verbose(0) << "Enter file name: ";
    1545           cin >> filename;
    1546         } while (!mol->AddXYZFile(filename));
    1547         break;
    1548 
    1549       case 'q': // quit
    1550         break;
    1551 
    1552       case 'r': // remove atom
    1553         RemoveAtoms(mol);
    1554         break;
    1555 
    1556       case 's': // save to config file
    1557         SaveConfig(ConfigFileName, &configuration, periode, mol);
    1558         break;
    1559 
    1560       case 't': // translate all atoms
    1561       cout << Verbose(0) << "Enter translation vector." << endl;
    1562       x.AskPosition(mol->cell_size,0);
    1563       mol->Translate((const Vector *)&x);
    1564       break;
    1565 
    1566       case 'T':
    1567         testroutine(mol);
    1568         break;
    1569 
    1570       case 'u': // change an atom's element
    1571         first = NULL;
    1572         do {
    1573           cout << Verbose(0) << "Change the element of which atom: ";
    1574           cin >> Z;
    1575         } while ((first = mol->FindAtom(Z)) == NULL);
    1576         cout << Verbose(0) << "New element by atomic number Z: ";
    1577         cin >> Z;
    1578         first->type = periode->FindElement(Z);
    1579         cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
    1580         break;
    1581     };
    1582   } while (choice != 'q');
    1583 
    1584   // save element data base
    1585   if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
    1586     cout << Verbose(0) << "Saving of elements.db successful." << endl;
    1587   else
    1588     cout << Verbose(0) << "Saving of elements.db failed." << endl;
    1589 
    1590   // Free all
    1591   if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
    1592     while (Subgraphs->next != NULL) {
    1593       Subgraphs = Subgraphs->next;
    1594       delete(Subgraphs->previous);
    1595     }
    1596     delete(Subgraphs);
    1597   }
    1598   delete(mol);
    1599   delete(periode);
    1600   return (0);
     1347        bool valid; // flag if input was valid or not
     1348        ifstream test;
     1349        ofstream output;
     1350        string line;
     1351        char filename[MAXSTRINGSIZE];
     1352        char *ConfigFileName = NULL;
     1353        char *ElementsFileName = NULL;
     1354        int Z;
     1355        int j, axis, count, faktor;
     1356        clock_t start,end;
     1357//      int *MinimumRingSize = NULL;
     1358        MoleculeLeafClass *Subgraphs = NULL;
     1359//      class StackClass<bond *> *BackEdgeStack = NULL;
     1360        element **Elements;
     1361        Vector **vectors;
     1362
     1363        // =========================== PARSE COMMAND LINE OPTIONS ====================================
     1364        j = ParseCommandLineOptions(argc, argv, mol, periode, configuration, ConfigFileName, ElementsFileName);
     1365        if (j == 1) return 0; // just for -v and -h options
     1366        if (j) return j;        // something went wrong
     1367
     1368        // General stuff
     1369        if (mol->cell_size[0] == 0.) {
     1370                cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
     1371                for (int i=0;i<6;i++) {
     1372                        cout << Verbose(1) << "Cell size" << i << ": ";
     1373                        cin >> mol->cell_size[i];
     1374                }
     1375        }
     1376
     1377        // =========================== START INTERACTIVE SESSION ====================================
     1378
     1379        // now the main construction loop
     1380        cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
     1381        do {
     1382                cout << Verbose(0) << endl << endl;
     1383                cout << Verbose(0) << "============Element list=======================" << endl;
     1384                mol->Checkout((ofstream *)&cout);
     1385                cout << Verbose(0) << "============Atom list==========================" << endl;
     1386                mol->Output((ofstream *)&cout);
     1387                cout << Verbose(0) << "============Menu===============================" << endl;
     1388                cout << Verbose(0) << "a - add an atom" << endl;
     1389                cout << Verbose(0) << "r - remove an atom" << endl;
     1390                cout << Verbose(0) << "b - scale a bond between atoms" << endl;
     1391                cout << Verbose(0) << "u - change an atoms element" << endl;
     1392                cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
     1393                cout << Verbose(0) << "-----------------------------------------------" << endl;
     1394                cout << Verbose(0) << "p - Parse xyz file" << endl;
     1395                cout << Verbose(0) << "e - edit the current configuration" << endl;
     1396                cout << Verbose(0) << "o - create connection matrix" << endl;
     1397                cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
     1398                cout << Verbose(0) << "-----------------------------------------------" << endl;
     1399                cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
     1400                cout << Verbose(0) << "i - realign molecule" << endl;
     1401                cout << Verbose(0) << "m - mirror all molecules" << endl;
     1402                cout << Verbose(0) << "t - translate molecule by vector" << endl;
     1403                cout << Verbose(0) << "c - scale by unit transformation" << endl;
     1404                cout << Verbose(0) << "g - center atoms in box" << endl;
     1405                cout << Verbose(0) << "-----------------------------------------------" << endl;
     1406                cout << Verbose(0) << "s - save current setup to config file" << endl;
     1407                cout << Verbose(0) << "T - call the current test routine" << endl;
     1408                cout << Verbose(0) << "q - quit" << endl;
     1409                cout << Verbose(0) << "===============================================" << endl;
     1410                cout << Verbose(0) << "Input: ";
     1411                cin >> choice;
     1412
     1413                switch (choice) {
     1414                        default:
     1415                        case 'a': // add atom
     1416                                AddAtoms(periode, mol);
     1417                                choice = 'a';
     1418                                break;
     1419
     1420                        case 'b': // scale a bond
     1421                                cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
     1422                                first = mol->AskAtom("Enter first (fixed) atom: ");
     1423                                second = mol->AskAtom("Enter second (shifting) atom: ");
     1424                                min_bond = 0.;
     1425                                for (int i=NDIM;i--;)
     1426                                        min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
     1427                                min_bond = sqrt(min_bond);
     1428                                cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
     1429                                cout << Verbose(0) << "Enter new bond length [a.u.]: ";
     1430                                cin >> bond;
     1431                                for (int i=NDIM;i--;) {
     1432                                        second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
     1433                                }
     1434                                //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
     1435                                //second->Output(second->type->No, 1, (ofstream *)&cout);
     1436                                break;
     1437
     1438                        case 'c': // unit scaling of the metric
     1439                        cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
     1440                        cout << Verbose(0) << "Enter three factors: ";
     1441                        factor = new double[NDIM];
     1442                        cin >> factor[0];
     1443                        cin >> factor[1];
     1444                        cin >> factor[2];
     1445                        valid = true;
     1446                        mol->Scale(&factor);
     1447                        delete[](factor);
     1448                        break;
     1449
     1450                        case 'd': // duplicate the periodic cell along a given axis, given times
     1451                                cout << Verbose(0) << "State the axis [(+-)123]: ";
     1452                                cin >> axis;
     1453                                cout << Verbose(0) << "State the factor: ";
     1454                                cin >> faktor;
     1455
     1456                                mol->CountAtoms((ofstream *)&cout);     // recount atoms
     1457                                if (mol->AtomCount != 0) {      // if there is more than none
     1458                                        count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
     1459                                        Elements = new element *[count];
     1460                                        vectors = new Vector *[count];
     1461                                        j = 0;
     1462                                        first = mol->start;
     1463                                        while (first->next != mol->end) {       // make a list of all atoms with coordinates and element
     1464                                                first = first->next;
     1465                                                Elements[j] = first->type;
     1466                                                vectors[j] = &first->x;
     1467                                                j++;
     1468                                        }
     1469                                        if (count != j)
     1470                                                cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
     1471                                        x.Zero();
     1472                                        y.Zero();
     1473                                        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
     1474                                        for (int i=1;i<faktor;i++) {    // then add this list with respective translation factor times
     1475                                                x.AddVector(&y); // per factor one cell width further
     1476                                                for (int k=count;k--;) { // go through every atom of the original cell
     1477                                                        first = new atom(); // create a new body
     1478                                                        first->x.CopyVector(vectors[k]);        // use coordinate of original atom
     1479                                                        first->x.AddVector(&x);                 // translate the coordinates
     1480                                                        first->type = Elements[k];      // insert original element
     1481                                                        mol->AddAtom(first);                            // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
     1482                                                }
     1483                                        }
     1484                                        if (mol->first->next != mol->last) // if connect matrix is present already, redo it
     1485                                                mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration.GetIsAngstroem());
     1486                                        // free memory
     1487                                        delete[](Elements);
     1488                                        delete[](vectors);
     1489                                        // correct cell size
     1490                                        if (axis < 0) { // if sign was negative, we have to translate everything
     1491                                                x.Zero();
     1492                                                x.AddVector(&y);
     1493                                                x.Scale(-(faktor-1));
     1494                                                mol->Translate(&x);
     1495                                        }
     1496                                        mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     1497                                }
     1498                                break;
     1499
     1500                        case 'e': // edit each field of the configuration
     1501                        configuration.Edit(mol);
     1502                        break;
     1503
     1504                        case 'f':
     1505                                FragmentAtoms(mol, &configuration);
     1506                                break;
     1507
     1508                        case 'g': // center the atoms
     1509                                CenterAtoms(mol);
     1510                                break;
     1511
     1512                        case 'i': // align all atoms
     1513                                AlignAtoms(periode, mol);
     1514                                break;
     1515
     1516                        case 'l': // measure distances or angles
     1517                                MeasureAtoms(periode, mol, &configuration);
     1518                                break;
     1519
     1520                        case 'm': // mirror atoms along a given axis
     1521                                MirrorAtoms(mol);
     1522                                break;
     1523
     1524                        case 'o': // create the connection matrix
     1525                                {
     1526                                        cout << Verbose(0) << "What's the maximum bond distance: ";
     1527                                        cin >> tmp1;
     1528                                        start = clock();
     1529                                        mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
     1530                                        mol->CreateListOfBondsPerAtom((ofstream *)&cout);
     1531//                                      Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
     1532//                                      while (Subgraphs->next != NULL) {
     1533//                                              Subgraphs = Subgraphs->next;
     1534//                                              Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
     1535//                                              delete(Subgraphs->previous);
     1536//                                      }
     1537//                                      delete(Subgraphs);              // we don't need the list here, so free everything
     1538//                                      delete[](MinimumRingSize);
     1539//                                      Subgraphs = NULL;
     1540                                        end = clock();
     1541                                        cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     1542                                }
     1543                                break;
     1544
     1545                        case 'p': // parse and XYZ file
     1546                                cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     1547                                do {
     1548                                        cout << Verbose(0) << "Enter file name: ";
     1549                                        cin >> filename;
     1550                                } while (!mol->AddXYZFile(filename));
     1551                                break;
     1552
     1553                        case 'q': // quit
     1554                                break;
     1555
     1556                        case 'r': // remove atom
     1557                                RemoveAtoms(mol);
     1558                                break;
     1559
     1560                        case 's': // save to config file
     1561                                SaveConfig(ConfigFileName, &configuration, periode, mol);
     1562                                break;
     1563
     1564                        case 't': // translate all atoms
     1565                        cout << Verbose(0) << "Enter translation vector." << endl;
     1566                        x.AskPosition(mol->cell_size,0);
     1567                        mol->Translate((const Vector *)&x);
     1568                        break;
     1569
     1570                        case 'T':
     1571                                testroutine(mol);
     1572                                break;
     1573
     1574                        case 'u': // change an atom's element
     1575                                first = NULL;
     1576                                do {
     1577                                        cout << Verbose(0) << "Change the element of which atom: ";
     1578                                        cin >> Z;
     1579                                } while ((first = mol->FindAtom(Z)) == NULL);
     1580                                cout << Verbose(0) << "New element by atomic number Z: ";
     1581                                cin >> Z;
     1582                                first->type = periode->FindElement(Z);
     1583                                cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
     1584                                break;
     1585                };
     1586        } while (choice != 'q');
     1587
     1588        // save element data base
     1589        if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
     1590                cout << Verbose(0) << "Saving of elements.db successful." << endl;
     1591        else
     1592                cout << Verbose(0) << "Saving of elements.db failed." << endl;
     1593
     1594        // Free all
     1595        if (Subgraphs != NULL) {        // free disconnected subgraph list of DFS analysis was performed
     1596                while (Subgraphs->next != NULL) {
     1597                        Subgraphs = Subgraphs->next;
     1598                        delete(Subgraphs->previous);
     1599                }
     1600                delete(Subgraphs);
     1601        }
     1602        delete(mol);
     1603        delete(periode);
     1604        return (0);
    16011605}
    16021606
Note: See TracChangeset for help on using the changeset viewer.