Changeset 6ac7ee for src/config.cpp


Ignore:
Timestamp:
Feb 9, 2009, 5:24:10 PM (17 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, Candidate_v1.7.0, 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/config.cpp

    • Property mode changed from 100644 to 100755
    r124df1 r6ac7ee  
    1313config::config()
    1414{
    15   mainname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    16   defaultpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    17   pseudopotpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    18   configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    19   configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    20   strcpy(mainname,"pcp");
    21   strcpy(defaultpath,"not specified");
    22   strcpy(pseudopotpath,"not specified");
    23   configpath[0]='\0';
    24   configname[0]='\0';
    25  
    26   FastParsing = false;
    27   ProcPEGamma=8;
    28   ProcPEPsi=1;
    29   DoOutVis=0;
    30   DoOutMes=1;
    31   DoOutNICS=0;
    32   DoOutOrbitals=0;
    33   DoOutCurrent=0;
    34   DoPerturbation=0;
    35   DoFullCurrent=0;
    36   DoWannier=0;
    37   CommonWannier=0;
    38   SawtoothStart=0.01;
    39   VectorPlane=0;
    40   VectorCut=0;
    41   UseAddGramSch=1;
    42   Seed=1;
    43  
    44   MaxOuterStep=0;
    45   Deltat=1;
    46   OutVisStep=10;
    47   OutSrcStep=5;
    48   TargetTemp=0.00095004455;
    49   ScaleTempStep=25;
    50   MaxPsiStep=0;
    51   EpsWannier=1e-7;
    52  
    53   MaxMinStep=100;
    54   RelEpsTotalEnergy=1e-7;
    55   RelEpsKineticEnergy=1e-5;
    56   MaxMinStopStep=1;
    57   MaxMinGapStopStep=0;
    58   MaxInitMinStep=100;
    59   InitRelEpsTotalEnergy=1e-5;
    60   InitRelEpsKineticEnergy=1e-4;
    61   InitMaxMinStopStep=1;
    62   InitMaxMinGapStopStep=0;
    63  
    64   //BoxLength[NDIM*NDIM];
    65  
    66   ECut=128.;
    67   MaxLevel=5;
    68   RiemannTensor=0;
    69   LevRFactor=0;
    70   RiemannLevel=0;
    71   Lev0Factor=2;
    72   RTActualUse=0;
    73   PsiType=0;
    74   MaxPsiDouble=0;
    75   PsiMaxNoUp=0;
    76   PsiMaxNoDown=0;
    77   AddPsis=0;
    78  
    79   RCut=20.;
    80   StructOpt=0;
    81   IsAngstroem=1;
    82   RelativeCoord=0;
    83   MaxTypes=0;
     15        mainname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     16        defaultpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     17        pseudopotpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     18        configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     19        configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     20        strcpy(mainname,"pcp");
     21        strcpy(defaultpath,"not specified");
     22        strcpy(pseudopotpath,"not specified");
     23        configpath[0]='\0';
     24        configname[0]='\0';
     25       
     26        FastParsing = false;
     27        ProcPEGamma=8;
     28        ProcPEPsi=1;
     29        DoOutVis=0;
     30        DoOutMes=1;
     31        DoOutNICS=0;
     32        DoOutOrbitals=0;
     33        DoOutCurrent=0;
     34        DoPerturbation=0;
     35        DoFullCurrent=0;
     36        DoWannier=0;
     37        CommonWannier=0;
     38        SawtoothStart=0.01;
     39        VectorPlane=0;
     40        VectorCut=0;
     41        UseAddGramSch=1;
     42        Seed=1;
     43       
     44        MaxOuterStep=0;
     45        Deltat=1;
     46        OutVisStep=10;
     47        OutSrcStep=5;
     48        TargetTemp=0.00095004455;
     49        ScaleTempStep=25;
     50        MaxPsiStep=0;
     51        EpsWannier=1e-7;
     52       
     53        MaxMinStep=100;
     54        RelEpsTotalEnergy=1e-7;
     55        RelEpsKineticEnergy=1e-5;
     56        MaxMinStopStep=1;
     57        MaxMinGapStopStep=0;
     58        MaxInitMinStep=100;
     59        InitRelEpsTotalEnergy=1e-5;
     60        InitRelEpsKineticEnergy=1e-4;
     61        InitMaxMinStopStep=1;
     62        InitMaxMinGapStopStep=0;
     63       
     64        //BoxLength[NDIM*NDIM];
     65       
     66        ECut=128.;
     67        MaxLevel=5;
     68        RiemannTensor=0;
     69        LevRFactor=0;
     70        RiemannLevel=0;
     71        Lev0Factor=2;
     72        RTActualUse=0;
     73        PsiType=0;
     74        MaxPsiDouble=0;
     75        PsiMaxNoUp=0;
     76        PsiMaxNoDown=0;
     77        AddPsis=0;
     78       
     79        RCut=20.;
     80        StructOpt=0;
     81        IsAngstroem=1;
     82        RelativeCoord=0;
     83        MaxTypes=0;
    8484};
    8585
     
    8989config::~config()
    9090{
    91   Free((void **)&mainname, "config::~config: *mainname");
    92   Free((void **)&defaultpath, "config::~config: *defaultpath");
    93   Free((void **)&pseudopotpath, "config::~config: *pseudopotpath");
    94   Free((void **)&configpath, "config::~config: *configpath");
    95   Free((void **)&configname, "config::~config: *configname");
     91        Free((void **)&mainname, "config::~config: *mainname");
     92        Free((void **)&defaultpath, "config::~config: *defaultpath");
     93        Free((void **)&pseudopotpath, "config::~config: *pseudopotpath");
     94        Free((void **)&configpath, "config::~config: *configpath");
     95        Free((void **)&configname, "config::~config: *configname");
    9696};
    9797
     
    102102void config::Edit(molecule *mol)
    103103{
    104   char choice;
    105  
    106   do {
    107     cout << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
    108     cout << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl;
    109     cout << Verbose(0) << " B - Default path (for runtime files)" << endl;
    110     cout << Verbose(0) << " C - Path of pseudopotential files" << endl;
    111     cout << Verbose(0) << " D - Number of coefficient sharing processes" << endl;
    112     cout << Verbose(0) << " E - Number of wave function sharing processes" << endl;
    113     cout << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl;
    114     cout << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl;
    115     cout << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl;
    116     cout << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl;
    117     cout << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl;
    118     cout << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl;
    119     cout << Verbose(0) << " L - 0: Wannier centres as calculated, 1: common centre for all, 2: unite centres according to spread, 3: cell centre, 4: shifted to nearest grid point" << endl;
    120     cout << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl;
    121     cout << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl;
    122     cout << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl;
    123     cout << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl;
    124     cout << Verbose(0) << " Q - Initial integer value of random number generator" << endl;
    125     cout << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl;
    126     cout << Verbose(0) << " T - Output visual after ...th step" << endl;
    127     cout << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl;
    128     cout << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl;
    129     cout << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl;
    130     cout << Verbose(0) << " Z - Maximum number of minimization iterations" << endl;
    131     cout << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl;
    132     cout << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl;
    133     cout << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl;
    134     cout << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl;
    135     cout << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl;
    136     cout << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
    137     cout << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
    138     cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
    139     cout << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
    140     cout << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
    141     cout << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl;
    142     cout << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl;
    143     cout << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl;
    144     cout << Verbose(0) << " p - Number of Riemann levels" << endl;
    145     cout << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl;
    146     cout << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl;
    147     cout << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl;
    148     cout << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl;
    149     cout << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl;
    150     cout << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl;
    151     cout << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl;
    152     cout << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl;
    153     cout << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl;
    154     cout << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl;
    155     cout << Verbose(0) << "=========================================================" << endl;
    156     cout << Verbose(0) << "INPUT: ";
    157     cin >> choice;
    158    
    159     switch (choice) {
    160         case 'A': // mainname
    161           cout << Verbose(0) << "Old: " << config::mainname << "\t new: ";
    162           cin >> config::mainname;
    163           break;
    164         case 'B': // defaultpath
    165           cout << Verbose(0) << "Old: " << config::defaultpath << "\t new: ";
    166           cin >> config::defaultpath;
    167           break;
    168         case 'C': // pseudopotpath
    169           cout << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ";
    170           cin >> config::pseudopotpath;
    171           break;
    172 
    173         case 'D': // ProcPEGamma
    174           cout << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ";
    175           cin >> config::ProcPEGamma;
    176           break;
    177         case 'E': // ProcPEPsi
    178           cout << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ";
    179           cin >> config::ProcPEPsi;
    180           break;
    181         case 'F': // DoOutVis
    182           cout << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ";
    183           cin >> config::DoOutVis;
    184           break;
    185         case 'G': // DoOutMes
    186           cout << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ";
    187           cin >> config::DoOutMes;
    188           break;
    189         case 'H': // DoOutOrbitals
    190           cout << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ";
    191           cin >> config::DoOutOrbitals;
    192           break;
    193         case 'I': // DoOutCurrent
    194           cout << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ";
    195           cin >> config::DoOutCurrent;
    196           break;
    197         case 'J': // DoFullCurrent
    198           cout << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ";
    199           cin >> config::DoFullCurrent;
    200           break;
    201         case 'K': // DoPerturbation
    202           cout << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ";
    203           cin >> config::DoPerturbation;
    204           break;
    205         case 'L': // CommonWannier
    206           cout << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ";
    207           cin >> config::CommonWannier;
    208           break;
    209         case 'M': // SawtoothStart
    210           cout << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ";
    211           cin >> config::SawtoothStart;
    212           break;
    213         case 'N': // VectorPlane
    214           cout << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ";
    215           cin >> config::VectorPlane;
    216           break;
    217         case 'O': // VectorCut
    218           cout << Verbose(0) << "Old: " << config::VectorCut << "\t new: ";
    219           cin >> config::VectorCut;
    220           break;
    221         case 'P': // UseAddGramSch
    222           cout << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ";
    223           cin >> config::UseAddGramSch;
    224           break;
    225         case 'Q': // Seed
    226           cout << Verbose(0) << "Old: " << config::Seed << "\t new: ";
    227           cin >> config::Seed;
    228           break;
    229 
    230         case 'R': // MaxOuterStep
    231           cout << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ";
    232           cin >> config::MaxOuterStep;
    233           break;
    234         case 'T': // OutVisStep
    235           cout << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ";
    236           cin >> config::OutVisStep;
    237           break;
    238         case 'U': // OutSrcStep
    239           cout << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ";
    240           cin >> config::OutSrcStep;
    241           break;
    242         case 'X': // MaxPsiStep
    243           cout << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ";
    244           cin >> config::MaxPsiStep;
    245           break;
    246         case 'Y': // EpsWannier
    247           cout << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ";
    248           cin >> config::EpsWannier;
    249           break;
    250 
    251         case 'Z': // MaxMinStep
    252           cout << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ";
    253           cin >> config::MaxMinStep;
    254           break;
    255         case 'a': // RelEpsTotalEnergy
    256           cout << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ";
    257           cin >> config::RelEpsTotalEnergy;
    258           break;
    259         case 'b': // RelEpsKineticEnergy
    260           cout << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ";
    261           cin >> config::RelEpsKineticEnergy;
    262           break;
    263         case 'c': // MaxMinStopStep
    264           cout << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ";
    265           cin >> config::MaxMinStopStep;
    266           break;
    267         case 'e': // MaxInitMinStep
    268           cout << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ";
    269           cin >> config::MaxInitMinStep;
    270           break;
    271         case 'f': // InitRelEpsTotalEnergy
    272           cout << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ";
    273           cin >> config::InitRelEpsTotalEnergy;
    274           break;
    275         case 'g': // InitRelEpsKineticEnergy
    276           cout << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ";
    277           cin >> config::InitRelEpsKineticEnergy;
    278           break;
    279         case 'h': // InitMaxMinStopStep
    280           cout << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ";
    281           cin >> config::InitMaxMinStopStep;
    282           break;
    283        
    284         case 'j': // BoxLength
    285           cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
    286           for (int i=0;i<6;i++) {
    287             cout << Verbose(0) << "Cell size" << i << ": ";
    288             cin >> mol->cell_size[i];
    289           }
    290           break;
    291        
    292         case 'k': // ECut
    293           cout << Verbose(0) << "Old: " << config::ECut << "\t new: ";
    294           cin >> config::ECut;
    295           break;
    296         case 'l': // MaxLevel
    297           cout << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ";
    298           cin >> config::MaxLevel;
    299           break;
    300         case 'm': // RiemannTensor
    301           cout << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ";
    302           cin >> config::RiemannTensor;
    303           break;
    304         case 'n': // LevRFactor
    305           cout << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ";
    306           cin >> config::LevRFactor;
    307           break;
    308         case 'o': // RiemannLevel
    309           cout << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ";
    310           cin >> config::RiemannLevel;
    311           break;
    312         case 'p': // Lev0Factor
    313           cout << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ";
    314           cin >> config::Lev0Factor;
    315           break;
    316         case 'r': // RTActualUse
    317           cout << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ";
    318           cin >> config::RTActualUse;
    319           break;
    320         case 's': // PsiType
    321           cout << Verbose(0) << "Old: " << config::PsiType << "\t new: ";
    322           cin >> config::PsiType;
    323           break;
    324         case 't': // MaxPsiDouble
    325           cout << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ";
    326           cin >> config::MaxPsiDouble;
    327           break;
    328         case 'u': // PsiMaxNoUp
    329           cout << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ";
    330           cin >> config::PsiMaxNoUp;
    331           break;
    332         case 'v': // PsiMaxNoDown
    333           cout << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ";
    334           cin >> config::PsiMaxNoDown;
    335           break;
    336         case 'w': // AddPsis
    337           cout << Verbose(0) << "Old: " << config::AddPsis << "\t new: ";
    338           cin >> config::AddPsis;
    339           break;
    340 
    341         case 'x': // RCut
    342           cout << Verbose(0) << "Old: " << config::RCut << "\t new: ";
    343           cin >> config::RCut;
    344           break;
    345         case 'y': // StructOpt
    346           cout << Verbose(0) << "Old: " << config::StructOpt << "\t new: ";
    347           cin >> config::StructOpt;
    348           break;
    349         case 'z': // IsAngstroem
    350           cout << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ";
    351           cin >> config::IsAngstroem;
    352           break;
    353         case 'i': // RelativeCoord
    354           cout << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ";
    355           cin >> config::RelativeCoord;
    356           break;
    357     };
    358   } while (choice != 'q');
     104        char choice;
     105       
     106        do {
     107                cout << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
     108                cout << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl;
     109                cout << Verbose(0) << " B - Default path (for runtime files)" << endl;
     110                cout << Verbose(0) << " C - Path of pseudopotential files" << endl;
     111                cout << Verbose(0) << " D - Number of coefficient sharing processes" << endl;
     112                cout << Verbose(0) << " E - Number of wave function sharing processes" << endl;
     113                cout << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl;
     114                cout << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl;
     115                cout << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl;
     116                cout << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl;
     117                cout << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl;
     118                cout << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl;
     119                cout << Verbose(0) << " L - 0: Wannier centres as calculated, 1: common centre for all, 2: unite centres according to spread, 3: cell centre, 4: shifted to nearest grid point" << endl;
     120                cout << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl;
     121                cout << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl;
     122                cout << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl;
     123                cout << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl;
     124                cout << Verbose(0) << " Q - Initial integer value of random number generator" << endl;
     125                cout << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl;
     126                cout << Verbose(0) << " T - Output visual after ...th step" << endl;
     127                cout << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl;
     128                cout << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl;
     129                cout << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl;
     130                cout << Verbose(0) << " Z - Maximum number of minimization iterations" << endl;
     131                cout << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl;
     132                cout << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl;
     133                cout << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl;
     134                cout << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl;
     135                cout << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl;
     136                cout << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
     137                cout << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
     138                cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
     139                cout << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
     140                cout << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
     141                cout << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl;
     142                cout << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl;
     143                cout << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl;
     144                cout << Verbose(0) << " p - Number of Riemann levels" << endl;
     145                cout << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl;
     146                cout << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl;
     147                cout << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl;
     148                cout << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl;
     149                cout << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl;
     150                cout << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl;
     151                cout << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl;
     152                cout << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl;
     153                cout << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl;
     154                cout << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl;
     155                cout << Verbose(0) << "=========================================================" << endl;
     156                cout << Verbose(0) << "INPUT: ";
     157                cin >> choice;
     158               
     159                switch (choice) {
     160                                case 'A': // mainname
     161                                        cout << Verbose(0) << "Old: " << config::mainname << "\t new: ";
     162                                        cin >> config::mainname;
     163                                        break;
     164                                case 'B': // defaultpath
     165                                        cout << Verbose(0) << "Old: " << config::defaultpath << "\t new: ";
     166                                        cin >> config::defaultpath;
     167                                        break;
     168                                case 'C': // pseudopotpath
     169                                        cout << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ";
     170                                        cin >> config::pseudopotpath;
     171                                        break;
     172
     173                                case 'D': // ProcPEGamma
     174                                        cout << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ";
     175                                        cin >> config::ProcPEGamma;
     176                                        break;
     177                                case 'E': // ProcPEPsi
     178                                        cout << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ";
     179                                        cin >> config::ProcPEPsi;
     180                                        break;
     181                                case 'F': // DoOutVis
     182                                        cout << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ";
     183                                        cin >> config::DoOutVis;
     184                                        break;
     185                                case 'G': // DoOutMes
     186                                        cout << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ";
     187                                        cin >> config::DoOutMes;
     188                                        break;
     189                                case 'H': // DoOutOrbitals
     190                                        cout << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ";
     191                                        cin >> config::DoOutOrbitals;
     192                                        break;
     193                                case 'I': // DoOutCurrent
     194                                        cout << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ";
     195                                        cin >> config::DoOutCurrent;
     196                                        break;
     197                                case 'J': // DoFullCurrent
     198                                        cout << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ";
     199                                        cin >> config::DoFullCurrent;
     200                                        break;
     201                                case 'K': // DoPerturbation
     202                                        cout << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ";
     203                                        cin >> config::DoPerturbation;
     204                                        break;
     205                                case 'L': // CommonWannier
     206                                        cout << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ";
     207                                        cin >> config::CommonWannier;
     208                                        break;
     209                                case 'M': // SawtoothStart
     210                                        cout << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ";
     211                                        cin >> config::SawtoothStart;
     212                                        break;
     213                                case 'N': // VectorPlane
     214                                        cout << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ";
     215                                        cin >> config::VectorPlane;
     216                                        break;
     217                                case 'O': // VectorCut
     218                                        cout << Verbose(0) << "Old: " << config::VectorCut << "\t new: ";
     219                                        cin >> config::VectorCut;
     220                                        break;
     221                                case 'P': // UseAddGramSch
     222                                        cout << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ";
     223                                        cin >> config::UseAddGramSch;
     224                                        break;
     225                                case 'Q': // Seed
     226                                        cout << Verbose(0) << "Old: " << config::Seed << "\t new: ";
     227                                        cin >> config::Seed;
     228                                        break;
     229
     230                                case 'R': // MaxOuterStep
     231                                        cout << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ";
     232                                        cin >> config::MaxOuterStep;
     233                                        break;
     234                                case 'T': // OutVisStep
     235                                        cout << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ";
     236                                        cin >> config::OutVisStep;
     237                                        break;
     238                                case 'U': // OutSrcStep
     239                                        cout << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ";
     240                                        cin >> config::OutSrcStep;
     241                                        break;
     242                                case 'X': // MaxPsiStep
     243                                        cout << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ";
     244                                        cin >> config::MaxPsiStep;
     245                                        break;
     246                                case 'Y': // EpsWannier
     247                                        cout << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ";
     248                                        cin >> config::EpsWannier;
     249                                        break;
     250
     251                                case 'Z': // MaxMinStep
     252                                        cout << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ";
     253                                        cin >> config::MaxMinStep;
     254                                        break;
     255                                case 'a': // RelEpsTotalEnergy
     256                                        cout << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ";
     257                                        cin >> config::RelEpsTotalEnergy;
     258                                        break;
     259                                case 'b': // RelEpsKineticEnergy
     260                                        cout << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ";
     261                                        cin >> config::RelEpsKineticEnergy;
     262                                        break;
     263                                case 'c': // MaxMinStopStep
     264                                        cout << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ";
     265                                        cin >> config::MaxMinStopStep;
     266                                        break;
     267                                case 'e': // MaxInitMinStep
     268                                        cout << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ";
     269                                        cin >> config::MaxInitMinStep;
     270                                        break;
     271                                case 'f': // InitRelEpsTotalEnergy
     272                                        cout << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ";
     273                                        cin >> config::InitRelEpsTotalEnergy;
     274                                        break;
     275                                case 'g': // InitRelEpsKineticEnergy
     276                                        cout << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ";
     277                                        cin >> config::InitRelEpsKineticEnergy;
     278                                        break;
     279                                case 'h': // InitMaxMinStopStep
     280                                        cout << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ";
     281                                        cin >> config::InitMaxMinStopStep;
     282                                        break;
     283                               
     284                                case 'j': // BoxLength
     285                                        cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
     286                                        for (int i=0;i<6;i++) {
     287                                                cout << Verbose(0) << "Cell size" << i << ": ";
     288                                                cin >> mol->cell_size[i];
     289                                        }
     290                                        break;
     291                               
     292                                case 'k': // ECut
     293                                        cout << Verbose(0) << "Old: " << config::ECut << "\t new: ";
     294                                        cin >> config::ECut;
     295                                        break;
     296                                case 'l': // MaxLevel
     297                                        cout << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ";
     298                                        cin >> config::MaxLevel;
     299                                        break;
     300                                case 'm': // RiemannTensor
     301                                        cout << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ";
     302                                        cin >> config::RiemannTensor;
     303                                        break;
     304                                case 'n': // LevRFactor
     305                                        cout << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ";
     306                                        cin >> config::LevRFactor;
     307                                        break;
     308                                case 'o': // RiemannLevel
     309                                        cout << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ";
     310                                        cin >> config::RiemannLevel;
     311                                        break;
     312                                case 'p': // Lev0Factor
     313                                        cout << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ";
     314                                        cin >> config::Lev0Factor;
     315                                        break;
     316                                case 'r': // RTActualUse
     317                                        cout << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ";
     318                                        cin >> config::RTActualUse;
     319                                        break;
     320                                case 's': // PsiType
     321                                        cout << Verbose(0) << "Old: " << config::PsiType << "\t new: ";
     322                                        cin >> config::PsiType;
     323                                        break;
     324                                case 't': // MaxPsiDouble
     325                                        cout << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ";
     326                                        cin >> config::MaxPsiDouble;
     327                                        break;
     328                                case 'u': // PsiMaxNoUp
     329                                        cout << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ";
     330                                        cin >> config::PsiMaxNoUp;
     331                                        break;
     332                                case 'v': // PsiMaxNoDown
     333                                        cout << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ";
     334                                        cin >> config::PsiMaxNoDown;
     335                                        break;
     336                                case 'w': // AddPsis
     337                                        cout << Verbose(0) << "Old: " << config::AddPsis << "\t new: ";
     338                                        cin >> config::AddPsis;
     339                                        break;
     340
     341                                case 'x': // RCut
     342                                        cout << Verbose(0) << "Old: " << config::RCut << "\t new: ";
     343                                        cin >> config::RCut;
     344                                        break;
     345                                case 'y': // StructOpt
     346                                        cout << Verbose(0) << "Old: " << config::StructOpt << "\t new: ";
     347                                        cin >> config::StructOpt;
     348                                        break;
     349                                case 'z': // IsAngstroem
     350                                        cout << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ";
     351                                        cin >> config::IsAngstroem;
     352                                        break;
     353                                case 'i': // RelativeCoord
     354                                        cout << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ";
     355                                        cin >> config::RelativeCoord;
     356                                        break;
     357                };
     358        } while (choice != 'q');
    359359};
    360360
     
    367367int config::TestSyntax(char *filename, periodentafel *periode, molecule *mol)
    368368{
    369   int test;
    370   ifstream file(filename);
    371  
    372   // search file for keyword: ProcPEGamma (new syntax)
    373   if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
    374     file.close();
    375     return 1;
    376   }
    377   // search file for keyword: ProcsGammaPsi (old syntax)
    378   if (ParseForParameter(1,&file,"ProcsGammaPsi", 0, 1, 1, int_type, &test, 1, optional)) {
    379     file.close();
    380     return 0;
    381   }
    382   file.close();
    383   return -1;
     369        int test;
     370        ifstream file(filename);
     371       
     372        // search file for keyword: ProcPEGamma (new syntax)
     373        if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
     374                file.close();
     375                return 1;
     376        }
     377        // search file for keyword: ProcsGammaPsi (old syntax)
     378        if (ParseForParameter(1,&file,"ProcsGammaPsi", 0, 1, 1, int_type, &test, 1, optional)) {
     379                file.close();
     380                return 0;
     381        }
     382        file.close();
     383        return -1;
    384384}
    385385
     
    389389bool config::GetIsAngstroem() const
    390390{
    391   return (IsAngstroem == 1);
     391        return (IsAngstroem == 1);
    392392};
    393393
     
    397397char * config::GetDefaultPath() const
    398398{
    399   return defaultpath;
     399        return defaultpath;
    400400};
    401401
     
    406406void config::SetDefaultPath(const char *path)
    407407{
    408   strcpy(defaultpath, path);
     408        strcpy(defaultpath, path);
    409409};
    410410
     
    414414void config::RetrieveConfigPathAndName(string filename)
    415415{
    416   char *ptr = NULL;
    417   char *buffer = new char[MAXSTRINGSIZE];
    418   strncpy(buffer, filename.c_str(), MAXSTRINGSIZE);
    419   int last = -1;
    420   for(last=MAXSTRINGSIZE;last--;) {
    421     if (buffer[last] == '/')
    422       break;
    423   }
    424   if (last == -1) { // no path in front, set to local directory.
    425     strcpy(configpath, "./");
    426     ptr = buffer;
    427   } else {
    428     strncpy(configpath, buffer, last+1);
    429     ptr = &buffer[last+1];
    430     if (last < 254)
    431       configpath[last+1]='\0';
    432   }
    433   strcpy(configname, ptr);
    434   cout << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl;
    435   delete[](buffer);
     416        char *ptr = NULL;
     417        char *buffer = new char[MAXSTRINGSIZE];
     418        strncpy(buffer, filename.c_str(), MAXSTRINGSIZE);
     419        int last = -1;
     420        for(last=MAXSTRINGSIZE;last--;) {
     421                if (buffer[last] == '/')
     422                        break;
     423        }
     424        if (last == -1) { // no path in front, set to local directory.
     425                strcpy(configpath, "./");
     426                ptr = buffer;
     427        } else {
     428                strncpy(configpath, buffer, last+1);
     429                ptr = &buffer[last+1];
     430                if (last < 254)
     431                        configpath[last+1]='\0';
     432        }
     433        strcpy(configname, ptr);
     434        cout << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl;
     435        delete[](buffer);
    436436};
    437437
     
    444444void config::Load(char *filename, periodentafel *periode, molecule *mol)
    445445{
    446   ifstream *file = new ifstream(filename);
    447   if (file == NULL) {
    448     cerr << "ERROR: config file " << filename << " missing!" << endl;
    449     return;
    450   }
    451   RetrieveConfigPathAndName(filename);
    452   // ParseParameters
    453  
    454   /* Oeffne Hauptparameterdatei */
    455   int di;
    456   double BoxLength[9];
    457   string zeile;
    458   string dummy;
    459   element *elementhash[MAX_ELEMENTS];
    460   char name[MAX_ELEMENTS];
    461   char keyword[MAX_ELEMENTS];
    462   int Z, No[MAX_ELEMENTS];
    463   int verbose = 0;
    464   double value[3];
    465  
    466   /* Namen einlesen */
    467 
    468   ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
    469   ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
    470   ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
    471   ParseForParameter(verbose,file,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
    472   ParseForParameter(verbose,file,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical);
    473 
    474   if (!ParseForParameter(verbose,file,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional))
    475     config::Seed = 1;
    476 
    477   if(!ParseForParameter(verbose,file,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) {
    478     config::DoOutOrbitals = 0;
    479   } else {
    480     if (config::DoOutOrbitals < 0) config::DoOutOrbitals = 0;
    481     if (config::DoOutOrbitals > 1) config::DoOutOrbitals = 1;
    482   }
    483   ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
    484   if (config::DoOutVis < 0) config::DoOutVis = 0;
    485   if (config::DoOutVis > 1) config::DoOutVis = 1;
    486   if (!ParseForParameter(verbose,file,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional))
    487     config::VectorPlane = -1;
    488   if (!ParseForParameter(verbose,file,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional))
    489     config::VectorCut = 0.;
    490   ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
    491   if (config::DoOutMes < 0) config::DoOutMes = 0;
    492   if (config::DoOutMes > 1) config::DoOutMes = 1;
    493   if (!ParseForParameter(verbose,file,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional))
    494     config::DoOutCurrent = 0;
    495   if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
    496   if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
    497   ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
    498   if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
    499   if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
    500   if(!ParseForParameter(verbose,file,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) {
    501     config::DoWannier = 0;
    502   } else {
    503     if (config::DoWannier < 0) config::DoWannier = 0;
    504     if (config::DoWannier > 1) config::DoWannier = 1;
    505   }
    506   if(!ParseForParameter(verbose,file,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) {
    507     config::CommonWannier = 0;
    508   } else {
    509     if (config::CommonWannier < 0) config::CommonWannier = 0;
    510     if (config::CommonWannier > 4) config::CommonWannier = 4;
    511   }
    512   if(!ParseForParameter(verbose,file,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) {
    513     config::SawtoothStart = 0.01;
    514   } else {
    515     if (config::SawtoothStart < 0.) config::SawtoothStart = 0.;
    516     if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
    517   }
    518  
    519   ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
    520   if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
    521     config::Deltat = 1;
    522   ParseForParameter(verbose,file,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    523   ParseForParameter(verbose,file,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    524   ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
    525   //ParseForParameter(verbose,file,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
    526   if (!ParseForParameter(verbose,file,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
    527     config::EpsWannier = 1e-8;
    528  
    529   // stop conditions
    530   //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
    531   ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
    532   if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
    533  
    534   ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
    535   ParseForParameter(verbose,file,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
    536   ParseForParameter(verbose,file,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
    537   ParseForParameter(verbose,file,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical);
    538   ParseForParameter(verbose,file,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical);
    539   if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
    540   if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
    541   if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;
    542  
    543   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
    544   ParseForParameter(verbose,file,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
    545   ParseForParameter(verbose,file,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
    546   ParseForParameter(verbose,file,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
    547   ParseForParameter(verbose,file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical);
    548   if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
    549   if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
    550   if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;
    551  
    552   // Unit cell and magnetic field
    553   ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    554   mol->cell_size[0] = BoxLength[0];
    555   mol->cell_size[1] = BoxLength[3];
    556   mol->cell_size[2] = BoxLength[4];
    557   mol->cell_size[3] = BoxLength[6];
    558   mol->cell_size[4] = BoxLength[7];
    559   mol->cell_size[5] = BoxLength[8];
    560   if (1) fprintf(stderr,"\n");
    561 
    562   ParseForParameter(verbose,file,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional);
    563   ParseForParameter(verbose,file,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional);
    564   if (!ParseForParameter(verbose,file,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional))
    565     config::DoFullCurrent = 0;
    566   if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
    567   if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
    568   if (config::DoOutNICS < 0) config::DoOutNICS = 0;
    569   if (config::DoOutNICS > 2) config::DoOutNICS = 2;
    570   if (config::DoPerturbation == 0) {
    571     config::DoFullCurrent = 0;
    572     config::DoOutNICS = 0;
    573   }
    574 
    575   ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
    576   ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
    577   ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
    578   if (config::Lev0Factor < 2) {
    579     config::Lev0Factor = 2;
    580   }
    581   ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
    582   if (di >= 0 && di < 2) {
    583     config::RiemannTensor = di;
    584   } else {
    585     fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    586     exit(1);
    587   }
    588   switch (config::RiemannTensor) {
    589     case 0: //UseNoRT
    590       if (config::MaxLevel < 2) {
    591         config::MaxLevel = 2;
    592       }
    593       config::LevRFactor = 2;
    594       config::RTActualUse = 0;
    595       break;
    596     case 1: // UseRT
    597       if (config::MaxLevel < 3) {
    598         config::MaxLevel = 3;
    599       }
    600       ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
    601       if (config::RiemannLevel < 2) {
    602         config::RiemannLevel = 2;
    603       }
    604       if (config::RiemannLevel > config::MaxLevel-1) {
    605         config::RiemannLevel = config::MaxLevel-1;
    606       }
    607       ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
    608       if (config::LevRFactor < 2) {
    609         config::LevRFactor = 2;
    610       }
    611       config::Lev0Factor = 2;
    612       config::RTActualUse = 2;
    613       break;
    614   }
    615   ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
    616   if (di >= 0 && di < 2) {
    617     config::PsiType = di;
    618   } else {
    619     fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    620     exit(1);
    621   }
    622   switch (config::PsiType) {
    623   case 0: // SpinDouble
    624     ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
    625     ParseForParameter(verbose,file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
    626     break;
    627   case 1: // SpinUpDown
    628     if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
    629     ParseForParameter(verbose,file,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
    630     ParseForParameter(verbose,file,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
    631     ParseForParameter(verbose,file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
    632     break;
    633   }
    634  
    635   // IonsInitRead
    636  
    637   ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
    638   ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
    639   ParseForParameter(verbose,file,"MaxTypes", 0, 1, 1, int_type, &(config::MaxTypes), 1, critical);
    640   if (!ParseForParameter(verbose,file,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional))
    641     config::RelativeCoord = 0;
    642   if (!ParseForParameter(verbose,file,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
    643     config::StructOpt = 0;
    644   if (MaxTypes == 0) {
    645     cerr << "There are no atoms according to MaxTypes in this config file." << endl;
    646   } else {
    647     // prescan number of ions per type
    648     cout << Verbose(0) << "Prescanning ions per type: " << endl;
    649     for (int i=0; i < config::MaxTypes; i++) {
    650       sprintf(name,"Ion_Type%i",i+1);
    651       ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
    652       ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
    653       elementhash[i] = periode->FindElement(Z);
    654       cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
    655     }
    656     int repetition = 0; // which repeated keyword shall be read
    657  
    658     map<int, atom *> AtomList[config::MaxTypes];
    659     if (!FastParsing) {
    660       // parse in trajectories
    661       bool status = true;
    662       atom *neues = NULL;
    663       while (status) {
    664         cout << "Currently parsing MD step " << repetition << "." << endl;
    665         for (int i=0; i < config::MaxTypes; i++) {
    666           sprintf(name,"Ion_Type%i",i+1);
    667           for(int j=0;j<No[i];j++) {
    668             sprintf(keyword,"%s_%i",name, j+1);
    669             if (repetition == 0) {
    670               neues = new atom();
    671               AtomList[i][j] = neues;
    672               neues->type = elementhash[i]; // find element type
    673               mol->AddAtom(neues);
    674             } else
    675               neues = AtomList[i][j];
    676             status = (status &&
    677                     ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
    678                     ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
    679                     ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
    680                     ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
    681             if (!status) break;
    682    
    683             // check size of vectors
    684             if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
    685               //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
    686               mol->Trajectories[neues].R.resize(repetition+10);
    687               mol->Trajectories[neues].U.resize(repetition+10);
    688               mol->Trajectories[neues].F.resize(repetition+10);
    689             }
    690          
    691             // put into trajectories list
    692             for (int d=0;d<NDIM;d++)
    693               mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
    694            
    695             // parse velocities if present
    696             if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
    697               neues->v.x[0] = 0.;
    698             if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
    699               neues->v.x[1] = 0.;
    700             if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
    701               neues->v.x[2] = 0.;
    702             for (int d=0;d<NDIM;d++)
    703               mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
    704      
    705             // parse forces if present
    706             if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
    707               value[0] = 0.;
    708             if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
    709               value[1] = 0.;
    710             if(!ParseForParameter(verbose,file, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
    711               value[2] = 0.;
    712             for (int d=0;d<NDIM;d++)
    713               mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
    714  
    715   //            cout << "Parsed position of step " << (repetition) << ": (";
    716   //            for (int d=0;d<NDIM;d++)
    717   //              cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";          // next step
    718   //            cout << ")\t(";
    719   //            for (int d=0;d<NDIM;d++)
    720   //              cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";          // next step
    721   //            cout << ")\t(";
    722   //            for (int d=0;d<NDIM;d++)
    723   //              cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";          // next step
    724   //            cout << ")" << endl;
    725           }
    726         }
    727         repetition++;
    728       }
    729       repetition--;
    730       cout << "Found " << repetition << " trajectory steps." << endl;
    731       mol->MDSteps = repetition;
    732     } else {
    733       // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)
    734       repetition = 0;
    735       while ( ParseForParameter(verbose,file, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
    736               ParseForParameter(verbose,file, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
    737               ParseForParameter(verbose,file, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
    738         repetition++;
    739       cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
    740       // parse in molecule coordinates
    741       for (int i=0; i < config::MaxTypes; i++) {
    742         sprintf(name,"Ion_Type%i",i+1);
    743         for(int j=0;j<No[i];j++) {
    744           sprintf(keyword,"%s_%i",name, j+1);
    745           atom *neues = new atom();
    746           // then parse for each atom the coordinates as often as present
    747           ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
    748           ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
    749           ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
    750           ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
    751           if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
    752             neues->v.x[0] = 0.;
    753           if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
    754             neues->v.x[1] = 0.;
    755           if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
    756             neues->v.x[2] = 0.;
    757           // here we don't care if forces are present (last in trajectories is always equal to current position)
    758           neues->type = elementhash[i]; // find element type
    759           mol->AddAtom(neues);
    760         }
    761       }
    762     }
    763   }
    764   file->close();
    765   delete(file);
     446        ifstream *file = new ifstream(filename);
     447        if (file == NULL) {
     448                cerr << "ERROR: config file " << filename << " missing!" << endl;
     449                return;
     450        }
     451        RetrieveConfigPathAndName(filename);
     452        // ParseParameters
     453       
     454        /* Oeffne Hauptparameterdatei */
     455        int di;
     456        double BoxLength[9];
     457        string zeile;
     458        string dummy;
     459        element *elementhash[MAX_ELEMENTS];
     460        char name[MAX_ELEMENTS];
     461        char keyword[MAX_ELEMENTS];
     462        int Z, No[MAX_ELEMENTS];
     463        int verbose = 0;
     464        double value[3];
     465       
     466        /* Namen einlesen */
     467
     468        ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
     469        ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
     470        ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
     471        ParseForParameter(verbose,file,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
     472        ParseForParameter(verbose,file,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical);
     473
     474        if (!ParseForParameter(verbose,file,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional))
     475                config::Seed = 1;
     476
     477        if(!ParseForParameter(verbose,file,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) {
     478                config::DoOutOrbitals = 0;
     479        } else {
     480                if (config::DoOutOrbitals < 0) config::DoOutOrbitals = 0;
     481                if (config::DoOutOrbitals > 1) config::DoOutOrbitals = 1;
     482        }
     483        ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
     484        if (config::DoOutVis < 0) config::DoOutVis = 0;
     485        if (config::DoOutVis > 1) config::DoOutVis = 1;
     486        if (!ParseForParameter(verbose,file,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional))
     487                config::VectorPlane = -1;
     488        if (!ParseForParameter(verbose,file,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional))
     489                config::VectorCut = 0.;
     490        ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
     491        if (config::DoOutMes < 0) config::DoOutMes = 0;
     492        if (config::DoOutMes > 1) config::DoOutMes = 1;
     493        if (!ParseForParameter(verbose,file,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional))
     494                config::DoOutCurrent = 0;
     495        if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
     496        if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
     497        ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
     498        if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
     499        if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
     500        if(!ParseForParameter(verbose,file,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) {
     501                config::DoWannier = 0;
     502        } else {
     503                if (config::DoWannier < 0) config::DoWannier = 0;
     504                if (config::DoWannier > 1) config::DoWannier = 1;
     505        }
     506        if(!ParseForParameter(verbose,file,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) {
     507                config::CommonWannier = 0;
     508        } else {
     509                if (config::CommonWannier < 0) config::CommonWannier = 0;
     510                if (config::CommonWannier > 4) config::CommonWannier = 4;
     511        }
     512        if(!ParseForParameter(verbose,file,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) {
     513                config::SawtoothStart = 0.01;
     514        } else {
     515                if (config::SawtoothStart < 0.) config::SawtoothStart = 0.;
     516                if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
     517        }
     518       
     519        ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
     520        if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
     521                config::Deltat = 1;
     522        ParseForParameter(verbose,file,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
     523        ParseForParameter(verbose,file,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
     524        ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
     525        //ParseForParameter(verbose,file,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
     526        if (!ParseForParameter(verbose,file,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
     527                config::EpsWannier = 1e-8;
     528       
     529        // stop conditions
     530        //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
     531        ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
     532        if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
     533       
     534        ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
     535        ParseForParameter(verbose,file,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
     536        ParseForParameter(verbose,file,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
     537        ParseForParameter(verbose,file,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical);
     538        ParseForParameter(verbose,file,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical);
     539        if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
     540        if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
     541        if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;
     542       
     543        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
     544        ParseForParameter(verbose,file,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
     545        ParseForParameter(verbose,file,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
     546        ParseForParameter(verbose,file,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
     547        ParseForParameter(verbose,file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical);
     548        if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
     549        if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
     550        if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;
     551       
     552        // Unit cell and magnetic field
     553        ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
     554        mol->cell_size[0] = BoxLength[0];
     555        mol->cell_size[1] = BoxLength[3];
     556        mol->cell_size[2] = BoxLength[4];
     557        mol->cell_size[3] = BoxLength[6];
     558        mol->cell_size[4] = BoxLength[7];
     559        mol->cell_size[5] = BoxLength[8];
     560        if (1) fprintf(stderr,"\n");
     561
     562        ParseForParameter(verbose,file,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional);
     563        ParseForParameter(verbose,file,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional);
     564        if (!ParseForParameter(verbose,file,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional))
     565                config::DoFullCurrent = 0;
     566        if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
     567        if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
     568        if (config::DoOutNICS < 0) config::DoOutNICS = 0;
     569        if (config::DoOutNICS > 2) config::DoOutNICS = 2;
     570        if (config::DoPerturbation == 0) {
     571                config::DoFullCurrent = 0;
     572                config::DoOutNICS = 0;
     573        }
     574
     575        ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
     576        ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
     577        ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
     578        if (config::Lev0Factor < 2) {
     579                config::Lev0Factor = 2;
     580        }
     581        ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
     582        if (di >= 0 && di < 2) {
     583                config::RiemannTensor = di;
     584        } else {
     585                fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
     586                exit(1);
     587        }
     588        switch (config::RiemannTensor) {
     589                case 0: //UseNoRT
     590                        if (config::MaxLevel < 2) {
     591                                config::MaxLevel = 2;
     592                        }
     593                        config::LevRFactor = 2;
     594                        config::RTActualUse = 0;
     595                        break;
     596                case 1: // UseRT
     597                        if (config::MaxLevel < 3) {
     598                                config::MaxLevel = 3;
     599                        }
     600                        ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
     601                        if (config::RiemannLevel < 2) {
     602                                config::RiemannLevel = 2;
     603                        }
     604                        if (config::RiemannLevel > config::MaxLevel-1) {
     605                                config::RiemannLevel = config::MaxLevel-1;
     606                        }
     607                        ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
     608                        if (config::LevRFactor < 2) {
     609                                config::LevRFactor = 2;
     610                        }
     611                        config::Lev0Factor = 2;
     612                        config::RTActualUse = 2;
     613                        break;
     614        }
     615        ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
     616        if (di >= 0 && di < 2) {
     617                config::PsiType = di;
     618        } else {
     619                fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
     620                exit(1);
     621        }
     622        switch (config::PsiType) {
     623        case 0: // SpinDouble
     624                ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
     625                ParseForParameter(verbose,file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
     626                break;
     627        case 1: // SpinUpDown
     628                if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
     629                ParseForParameter(verbose,file,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
     630                ParseForParameter(verbose,file,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
     631                ParseForParameter(verbose,file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
     632                break;
     633        }
     634       
     635        // IonsInitRead
     636       
     637        ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
     638        ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
     639        ParseForParameter(verbose,file,"MaxTypes", 0, 1, 1, int_type, &(config::MaxTypes), 1, critical);
     640        if (!ParseForParameter(verbose,file,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional))
     641                config::RelativeCoord = 0;
     642        if (!ParseForParameter(verbose,file,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
     643                config::StructOpt = 0;
     644        if (MaxTypes == 0) {
     645                cerr << "There are no atoms according to MaxTypes in this config file." << endl;
     646        } else {
     647                // prescan number of ions per type
     648                cout << Verbose(0) << "Prescanning ions per type: " << endl;
     649                for (int i=0; i < config::MaxTypes; i++) {
     650                        sprintf(name,"Ion_Type%i",i+1);
     651                        ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
     652                        ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
     653                        elementhash[i] = periode->FindElement(Z);
     654                        cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
     655                }
     656                int repetition = 0; // which repeated keyword shall be read
     657       
     658                map<int, atom *> AtomList[config::MaxTypes];
     659                if (!FastParsing) {
     660                        // parse in trajectories
     661                        bool status = true;
     662                        atom *neues = NULL;
     663                        while (status) {
     664                                cout << "Currently parsing MD step " << repetition << "." << endl;
     665                                for (int i=0; i < config::MaxTypes; i++) {
     666                                        sprintf(name,"Ion_Type%i",i+1);
     667                                        for(int j=0;j<No[i];j++) {
     668                                                sprintf(keyword,"%s_%i",name, j+1);
     669                                                if (repetition == 0) {
     670                                                        neues = new atom();
     671                                                        AtomList[i][j] = neues;
     672                                                        neues->type = elementhash[i]; // find element type
     673                                                        mol->AddAtom(neues);
     674                                                } else
     675                                                        neues = AtomList[i][j];
     676                                                status = (status &&
     677                                                                                ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
     678                                                                                ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
     679                                                                                ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
     680                                                                                ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
     681                                                if (!status) break;
     682               
     683                                                // check size of vectors
     684                                                if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
     685                                                        //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
     686                                                        mol->Trajectories[neues].R.resize(repetition+10);
     687                                                        mol->Trajectories[neues].U.resize(repetition+10);
     688                                                        mol->Trajectories[neues].F.resize(repetition+10);
     689                                                }
     690                                       
     691                                                // put into trajectories list
     692                                                for (int d=0;d<NDIM;d++)
     693                                                        mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
     694                                               
     695                                                // parse velocities if present
     696                                                if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
     697                                                        neues->v.x[0] = 0.;
     698                                                if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
     699                                                        neues->v.x[1] = 0.;
     700                                                if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
     701                                                        neues->v.x[2] = 0.;
     702                                                for (int d=0;d<NDIM;d++)
     703                                                        mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
     704                       
     705                                                // parse forces if present
     706                                                if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
     707                                                        value[0] = 0.;
     708                                                if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
     709                                                        value[1] = 0.;
     710                                                if(!ParseForParameter(verbose,file, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
     711                                                        value[2] = 0.;
     712                                                for (int d=0;d<NDIM;d++)
     713                                                        mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
     714       
     715        //                                              cout << "Parsed position of step " << (repetition) << ": (";
     716        //                                              for (int d=0;d<NDIM;d++)
     717        //                                                      cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";                                  // next step
     718        //                                              cout << ")\t(";
     719        //                                              for (int d=0;d<NDIM;d++)
     720        //                                                      cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";                                  // next step
     721        //                                              cout << ")\t(";
     722        //                                              for (int d=0;d<NDIM;d++)
     723        //                                                      cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";                                  // next step
     724        //                                              cout << ")" << endl;
     725                                        }
     726                                }
     727                                repetition++;
     728                        }
     729                        repetition--;
     730                        cout << "Found " << repetition << " trajectory steps." << endl;
     731                        mol->MDSteps = repetition;
     732                } else {
     733                        // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)
     734                        repetition = 0;
     735                        while ( ParseForParameter(verbose,file, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
     736                                                        ParseForParameter(verbose,file, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
     737                                                        ParseForParameter(verbose,file, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
     738                                repetition++;
     739                        cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
     740                        // parse in molecule coordinates
     741                        for (int i=0; i < config::MaxTypes; i++) {
     742                                sprintf(name,"Ion_Type%i",i+1);
     743                                for(int j=0;j<No[i];j++) {
     744                                        sprintf(keyword,"%s_%i",name, j+1);
     745                                        atom *neues = new atom();
     746                                        // then parse for each atom the coordinates as often as present
     747                                        ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
     748                                        ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
     749                                        ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
     750                                        ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
     751                                        if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
     752                                                neues->v.x[0] = 0.;
     753                                        if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
     754                                                neues->v.x[1] = 0.;
     755                                        if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
     756                                                neues->v.x[2] = 0.;
     757                                        // here we don't care if forces are present (last in trajectories is always equal to current position)
     758                                        neues->type = elementhash[i]; // find element type
     759                                        mol->AddAtom(neues);
     760                                }
     761                        }
     762                }
     763        }
     764        file->close();
     765        delete(file);
    766766};
    767767
     
    773773void config::LoadOld(char *filename, periodentafel *periode, molecule *mol)
    774774{
    775   ifstream *file = new ifstream(filename);
    776   if (file == NULL) {
    777     cerr << "ERROR: config file " << filename << " missing!" << endl;
    778     return;
    779   }
    780   RetrieveConfigPathAndName(filename);
    781   // ParseParameters
    782  
    783   /* Oeffne Hauptparameterdatei */
    784   int l, i, di;
    785   double a,b;
    786   double BoxLength[9];
    787   string zeile;
    788   string dummy;
    789   element *elementhash[128];
    790   int Z, No, AtomNo, found;
    791   int verbose = 0;
    792  
    793   /* Namen einlesen */
    794 
    795   ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
    796   ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
    797   ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
    798   ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
    799   ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 2, 1, int_type, &(config::ProcPEPsi), 1, critical);
    800   config::Seed = 1;
    801   config::DoOutOrbitals = 0;
    802   ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
    803   if (config::DoOutVis < 0) config::DoOutVis = 0;
    804   if (config::DoOutVis > 1) config::DoOutVis = 1;
    805   config::VectorPlane = -1;
    806   config::VectorCut = 0.;
    807   ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
    808   if (config::DoOutMes < 0) config::DoOutMes = 0;
    809   if (config::DoOutMes > 1) config::DoOutMes = 1;
    810   config::DoOutCurrent = 0;
    811   ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
    812   if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
    813   if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
    814   config::CommonWannier = 0;
    815   config::SawtoothStart = 0.01;
    816 
    817   ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, double_type, &(config::MaxOuterStep), 1, critical);
    818   ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
    819   ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    820   ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    821   ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
    822   ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
    823   config::EpsWannier = 1e-8;
    824  
    825   // stop conditions
    826   //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
    827   ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
    828   if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
    829  
    830   ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
    831   ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
    832   ParseForParameter(verbose,file,"MaxMinStep", 0, 3, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
    833   ParseForParameter(verbose,file,"MaxMinStep", 0, 4, 1, int_type, &(config::MaxMinStopStep), 1, critical);
    834   if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
    835   if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
    836   config::MaxMinGapStopStep = 1;
    837  
    838   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
    839   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
    840   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 3, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
    841   ParseForParameter(verbose,file,"MaxInitMinStep", 0, 4, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
    842   if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
    843   if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
    844   config::InitMaxMinGapStopStep = 1;
    845 
    846   ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    847   mol->cell_size[0] = BoxLength[0];
    848   mol->cell_size[1] = BoxLength[3];
    849   mol->cell_size[2] = BoxLength[4];
    850   mol->cell_size[3] = BoxLength[6];
    851   mol->cell_size[4] = BoxLength[7];
    852   mol->cell_size[5] = BoxLength[8];
    853   if (1) fprintf(stderr,"\n");
    854   config::DoPerturbation = 0;
    855   config::DoFullCurrent = 0;
    856 
    857   ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
    858   ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
    859   ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
    860   if (config::Lev0Factor < 2) {
    861     config::Lev0Factor = 2;
    862   }
    863   ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
    864   if (di >= 0 && di < 2) {
    865     config::RiemannTensor = di;
    866   } else {
    867     fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    868     exit(1);
    869   }
    870   switch (config::RiemannTensor) {
    871     case 0: //UseNoRT
    872       if (config::MaxLevel < 2) {
    873         config::MaxLevel = 2;
    874       }
    875       config::LevRFactor = 2;
    876       config::RTActualUse = 0;
    877       break;
    878     case 1: // UseRT
    879       if (config::MaxLevel < 3) {
    880         config::MaxLevel = 3;
    881       }
    882       ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
    883       if (config::RiemannLevel < 2) {
    884         config::RiemannLevel = 2;
    885       }
    886       if (config::RiemannLevel > config::MaxLevel-1) {
    887         config::RiemannLevel = config::MaxLevel-1;
    888       }
    889       ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
    890       if (config::LevRFactor < 2) {
    891         config::LevRFactor = 2;
    892       }
    893       config::Lev0Factor = 2;
    894       config::RTActualUse = 2;
    895       break;
    896   }
    897   ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
    898   if (di >= 0 && di < 2) {
    899     config::PsiType = di;
    900   } else {
    901     fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    902     exit(1);
    903   }
    904   switch (config::PsiType) {
    905   case 0: // SpinDouble
    906     ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
    907     config::AddPsis = 0;
    908     break;
    909   case 1: // SpinUpDown
    910     if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
    911     ParseForParameter(verbose,file,"MaxPsiUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
    912     ParseForParameter(verbose,file,"MaxPsiDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
    913     config::AddPsis = 0;
    914     break;
    915   }
    916  
    917   // IonsInitRead
    918  
    919   ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
    920   ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
    921   config::RelativeCoord = 0;
    922   config::StructOpt = 0;
    923 
    924   // Routine from builder.cpp
    925  
    926  
    927   for (i=MAX_ELEMENTS;i--;)
    928     elementhash[i] = NULL;
    929   cout << Verbose(0) << "Parsing Ions ..." << endl;
    930   No=0;
    931   found = 0;
    932   while (getline(*file,zeile,'\n')) {
    933     if (zeile.find("Ions_Data") == 0) {
    934       cout << Verbose(1) << "found Ions_Data...begin parsing" << endl;
    935       found ++;
    936     }
    937     if (found > 0) {
    938       if (zeile.find("Ions_Data") == 0)
    939         getline(*file,zeile,'\n'); // read next line and parse this one
    940       istringstream input(zeile);
    941       input >> AtomNo;  // number of atoms
    942       input >> Z;      // atomic number
    943       input >> a;
    944       input >> l;
    945       input >> l;
    946       input >> b;    // element mass
    947       elementhash[No] = periode->FindElement(Z);
    948       cout << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:"  << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl;
    949       for(i=0;i<AtomNo;i++) {
    950         if (!getline(*file,zeile,'\n')) {// parse on and on
    951           cout << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl;
    952           // return 1;
    953         } else {
    954           //cout << Verbose(2) << "Reading line: " << zeile << endl;
    955         }
    956         istringstream input2(zeile);
    957         atom *neues = new atom();
    958         input2 >> neues->x.x[0]; // x
    959         input2 >> neues->x.x[1]; // y
    960         input2 >> neues->x.x[2]; // z
    961         input2 >> l;
    962         neues->type = elementhash[No]; // find element type
    963         mol->AddAtom(neues);
    964       }
    965       No++;
    966     }
    967   } 
    968   file->close();
    969   delete(file);
     775        ifstream *file = new ifstream(filename);
     776        if (file == NULL) {
     777                cerr << "ERROR: config file " << filename << " missing!" << endl;
     778                return;
     779        }
     780        RetrieveConfigPathAndName(filename);
     781        // ParseParameters
     782       
     783        /* Oeffne Hauptparameterdatei */
     784        int l, i, di;
     785        double a,b;
     786        double BoxLength[9];
     787        string zeile;
     788        string dummy;
     789        element *elementhash[128];
     790        int Z, No, AtomNo, found;
     791        int verbose = 0;
     792       
     793        /* Namen einlesen */
     794
     795        ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
     796        ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
     797        ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
     798        ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
     799        ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 2, 1, int_type, &(config::ProcPEPsi), 1, critical);
     800        config::Seed = 1;
     801        config::DoOutOrbitals = 0;
     802        ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
     803        if (config::DoOutVis < 0) config::DoOutVis = 0;
     804        if (config::DoOutVis > 1) config::DoOutVis = 1;
     805        config::VectorPlane = -1;
     806        config::VectorCut = 0.;
     807        ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
     808        if (config::DoOutMes < 0) config::DoOutMes = 0;
     809        if (config::DoOutMes > 1) config::DoOutMes = 1;
     810        config::DoOutCurrent = 0;
     811        ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
     812        if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
     813        if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
     814        config::CommonWannier = 0;
     815        config::SawtoothStart = 0.01;
     816
     817        ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, double_type, &(config::MaxOuterStep), 1, critical);
     818        ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
     819        ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
     820        ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
     821        ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
     822        ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
     823        config::EpsWannier = 1e-8;
     824       
     825        // stop conditions
     826        //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
     827        ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
     828        if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
     829       
     830        ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
     831        ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
     832        ParseForParameter(verbose,file,"MaxMinStep", 0, 3, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
     833        ParseForParameter(verbose,file,"MaxMinStep", 0, 4, 1, int_type, &(config::MaxMinStopStep), 1, critical);
     834        if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
     835        if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
     836        config::MaxMinGapStopStep = 1;
     837       
     838        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
     839        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
     840        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 3, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
     841        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 4, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
     842        if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
     843        if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
     844        config::InitMaxMinGapStopStep = 1;
     845
     846        ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
     847        mol->cell_size[0] = BoxLength[0];
     848        mol->cell_size[1] = BoxLength[3];
     849        mol->cell_size[2] = BoxLength[4];
     850        mol->cell_size[3] = BoxLength[6];
     851        mol->cell_size[4] = BoxLength[7];
     852        mol->cell_size[5] = BoxLength[8];
     853        if (1) fprintf(stderr,"\n");
     854        config::DoPerturbation = 0;
     855        config::DoFullCurrent = 0;
     856
     857        ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
     858        ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
     859        ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
     860        if (config::Lev0Factor < 2) {
     861                config::Lev0Factor = 2;
     862        }
     863        ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
     864        if (di >= 0 && di < 2) {
     865                config::RiemannTensor = di;
     866        } else {
     867                fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
     868                exit(1);
     869        }
     870        switch (config::RiemannTensor) {
     871                case 0: //UseNoRT
     872                        if (config::MaxLevel < 2) {
     873                                config::MaxLevel = 2;
     874                        }
     875                        config::LevRFactor = 2;
     876                        config::RTActualUse = 0;
     877                        break;
     878                case 1: // UseRT
     879                        if (config::MaxLevel < 3) {
     880                                config::MaxLevel = 3;
     881                        }
     882                        ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
     883                        if (config::RiemannLevel < 2) {
     884                                config::RiemannLevel = 2;
     885                        }
     886                        if (config::RiemannLevel > config::MaxLevel-1) {
     887                                config::RiemannLevel = config::MaxLevel-1;
     888                        }
     889                        ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
     890                        if (config::LevRFactor < 2) {
     891                                config::LevRFactor = 2;
     892                        }
     893                        config::Lev0Factor = 2;
     894                        config::RTActualUse = 2;
     895                        break;
     896        }
     897        ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
     898        if (di >= 0 && di < 2) {
     899                config::PsiType = di;
     900        } else {
     901                fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
     902                exit(1);
     903        }
     904        switch (config::PsiType) {
     905        case 0: // SpinDouble
     906                ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
     907                config::AddPsis = 0;
     908                break;
     909        case 1: // SpinUpDown
     910                if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
     911                ParseForParameter(verbose,file,"MaxPsiUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
     912                ParseForParameter(verbose,file,"MaxPsiDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
     913                config::AddPsis = 0;
     914                break;
     915        }
     916       
     917        // IonsInitRead
     918       
     919        ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
     920        ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
     921        config::RelativeCoord = 0;
     922        config::StructOpt = 0;
     923
     924        // Routine from builder.cpp
     925       
     926       
     927        for (i=MAX_ELEMENTS;i--;)
     928                elementhash[i] = NULL;
     929        cout << Verbose(0) << "Parsing Ions ..." << endl;
     930        No=0;
     931        found = 0;
     932        while (getline(*file,zeile,'\n')) {
     933                if (zeile.find("Ions_Data") == 0) {
     934                        cout << Verbose(1) << "found Ions_Data...begin parsing" << endl;
     935                        found ++;
     936                }
     937                if (found > 0) {
     938                        if (zeile.find("Ions_Data") == 0)
     939                                getline(*file,zeile,'\n'); // read next line and parse this one
     940                        istringstream input(zeile);
     941                        input >> AtomNo;        // number of atoms
     942                        input >> Z;                      // atomic number
     943                        input >> a;
     944                        input >> l;
     945                        input >> l;
     946                        input >> b;              // element mass
     947                        elementhash[No] = periode->FindElement(Z);
     948                        cout << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:"     << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl;
     949                        for(i=0;i<AtomNo;i++) {
     950                                if (!getline(*file,zeile,'\n')) {// parse on and on
     951                                        cout << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl;
     952                                        // return 1;
     953                                } else {
     954                                        //cout << Verbose(2) << "Reading line: " << zeile << endl;
     955                                }
     956                                istringstream input2(zeile);
     957                                atom *neues = new atom();
     958                                input2 >> neues->x.x[0]; // x
     959                                input2 >> neues->x.x[1]; // y
     960                                input2 >> neues->x.x[2]; // z
     961                                input2 >> l;
     962                                neues->type = elementhash[No]; // find element type
     963                                mol->AddAtom(neues);
     964                        }
     965                        No++;
     966                }
     967        }       
     968        file->close();
     969        delete(file);
    970970};
    971971
     
    977977bool config::Save(const char *filename, periodentafel *periode, molecule *mol) const
    978978{
    979   bool result = true;
     979        bool result = true;
    980980        // bring MaxTypes up to date
    981981        mol->CountElements();
    982   ofstream *output = NULL;
    983   output = new ofstream(filename, ios::out);
    984   if (output != NULL) {
    985     *output << "# ParallelCarParinello - main configuration file - created with molecuilder" << endl;
    986     *output << endl;
    987     *output << "mainname\t" << config::mainname << "\t# programm name (for runtime files)" << endl;
    988     *output << "defaultpath\t" << config::defaultpath << "\t# where to put files during runtime" << endl;
    989     *output << "pseudopotpath\t" << config::pseudopotpath << "\t# where to find pseudopotentials" << endl;
    990     *output << endl;
    991     *output << "ProcPEGamma\t" << config::ProcPEGamma << "\t# for parallel computing: share constants" << endl;
    992     *output << "ProcPEPsi\t" << config::ProcPEPsi << "\t# for parallel computing: share wave functions" << endl;
    993     *output << "DoOutVis\t" << config::DoOutVis << "\t# Output data for OpenDX" << endl;
    994     *output << "DoOutMes\t" << config::DoOutMes << "\t# Output data for measurements" << endl;
    995     *output << "DoOutOrbitals\t" << config::DoOutOrbitals << "\t# Output all Orbitals" << endl;
    996     *output << "DoOutCurr\t" << config::DoOutCurrent << "\t# Ouput current density for OpenDx" << endl;
    997     *output << "DoOutNICS\t" << config::DoOutNICS << "\t# Output Nucleus independent current shieldings" << endl;
    998     *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
    999     *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
    1000     *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    1001     *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
    1002     *output << "VectorPlane\t" << config::VectorPlane << "\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot" << endl;
    1003     *output << "VectorCut\t" << config::VectorCut << "\t# Cut plane axis value" << endl;
    1004     *output << "AddGramSch\t" << config::UseAddGramSch << "\t# Additional GramSchmidtOrtogonalization to be safe" << endl;
    1005     *output << "Seed\t\t" << config::Seed << "\t# initial value for random seed for Psi coefficients" << endl;
    1006     *output << endl;
    1007     *output << "MaxOuterStep\t" << config::MaxOuterStep << "\t# number of MolecularDynamics/Structure optimization steps" << endl;
    1008     *output << "Deltat\t" << config::Deltat << "\t# time per MD step" << endl;
    1009     *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
    1010     *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
    1011     *output << "TargetTemp\t" << config::TargetTemp << "\t# Target temperature" << endl;
    1012     *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
    1013     *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
    1014     *output << endl;
    1015     *output << "# Values specifying when to stop" << endl;
    1016     *output << "MaxMinStep\t" << config::MaxMinStep << "\t# Maximum number of steps" << endl;
    1017     *output << "RelEpsTotalE\t" << config::RelEpsTotalEnergy << "\t# relative change in total energy" << endl;
    1018     *output << "RelEpsKineticE\t" << config::RelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
    1019     *output << "MaxMinStopStep\t" << config::MaxMinStopStep << "\t# check every ..th steps" << endl;
    1020     *output << "MaxMinGapStopStep\t" << config::MaxMinGapStopStep << "\t# check every ..th steps" << endl;
    1021     *output << endl;
    1022     *output << "# Values specifying when to stop for INIT, otherwise same as above" << endl;
    1023     *output << "MaxInitMinStep\t" << config::MaxInitMinStep << "\t# Maximum number of steps" << endl;
    1024     *output << "InitRelEpsTotalE\t" << config::InitRelEpsTotalEnergy << "\t# relative change in total energy" << endl;
    1025     *output << "InitRelEpsKineticE\t" << config::InitRelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
    1026     *output << "InitMaxMinStopStep\t" << config::InitMaxMinStopStep << "\t# check every ..th steps" << endl;
    1027     *output << "InitMaxMinGapStopStep\t" << config::InitMaxMinGapStopStep << "\t# check every ..th steps" << endl;
    1028     *output << endl;
    1029     *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    1030     *output << mol->cell_size[0] << "\t" << endl;
    1031     *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
    1032     *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
    1033     // FIXME
    1034     *output << endl;
    1035     *output << "ECut\t\t" << config::ECut << "\t# energy cutoff for discretization in Hartrees" << endl;
    1036     *output << "MaxLevel\t" << config::MaxLevel << "\t# number of different levels in the code, >=2" << endl;
    1037     *output << "Level0Factor\t" << config::Lev0Factor << "\t# factor by which node number increases from S to 0 level" << endl;
    1038     *output << "RiemannTensor\t" << config::RiemannTensor << "\t# (Use metric)" << endl;
    1039     switch (config::RiemannTensor) {
    1040       case 0: //UseNoRT
    1041         break;
    1042       case 1: // UseRT
    1043         *output << "RiemannLevel\t" << config::RiemannLevel << "\t# Number of Riemann Levels" << endl;
    1044         *output << "LevRFactor\t" << config::LevRFactor << "\t# factor by which node number increases from 0 to R level from" << endl;
    1045         break;
    1046     }
    1047     *output << "PsiType\t\t" << config::PsiType << "\t# 0 - doubly occupied, 1 - SpinUp,SpinDown" << endl;
    1048     // write out both types for easier changing afterwards
    1049   //  switch (PsiType) {
    1050   //    case 0:
    1051         *output << "MaxPsiDouble\t" << config::MaxPsiDouble << "\t# here: specifying both maximum number of SpinUp- and -Down-states" << endl;
    1052   //      break;
    1053   //    case 1:
    1054         *output << "PsiMaxNoUp\t" << config::PsiMaxNoUp << "\t# here: specifying maximum number of SpinUp-states" << endl;
    1055         *output << "PsiMaxNoDown\t" << config::PsiMaxNoDown << "\t# here: specifying maximum number of SpinDown-states" << endl;
    1056   //      break;
    1057   //  }
    1058     *output << "AddPsis\t\t" << config::AddPsis << "\t# Additional unoccupied Psis for bandgap determination" << endl;
    1059     *output << endl;
    1060     *output << "RCut\t\t" << config::RCut << "\t# R-cut for the ewald summation" << endl;
    1061     *output << "StructOpt\t" << config::StructOpt << "\t# Do structure optimization beforehand" << endl;
    1062     *output << "IsAngstroem\t" << config::IsAngstroem << "\t# 0 - Bohr, 1 - Angstroem" << endl;
    1063     *output << "RelativeCoord\t" << config::RelativeCoord << "\t# whether ion coordinates are relative (1) or absolute (0)" << endl;
    1064     *output << "MaxTypes\t" << mol->ElementCount <<  "\t# maximum number of different ion types" << endl;
    1065     *output << endl;
    1066     result = result && mol->Checkout(output);
    1067     if (mol->MDSteps <=1 )
    1068       result = result && mol->Output(output);
    1069     else
    1070       result = result && mol->OutputTrajectories(output);
    1071     output->close();
    1072     output->clear();
    1073     delete(output);
    1074     return result;
    1075   } else
    1076     return false;
     982        ofstream *output = NULL;
     983        output = new ofstream(filename, ios::out);
     984        if (output != NULL) {
     985                *output << "# ParallelCarParinello - main configuration file - created with molecuilder" << endl;
     986                *output << endl;
     987                *output << "mainname\t" << config::mainname << "\t# programm name (for runtime files)" << endl;
     988                *output << "defaultpath\t" << config::defaultpath << "\t# where to put files during runtime" << endl;
     989                *output << "pseudopotpath\t" << config::pseudopotpath << "\t# where to find pseudopotentials" << endl;
     990                *output << endl;
     991                *output << "ProcPEGamma\t" << config::ProcPEGamma << "\t# for parallel computing: share constants" << endl;
     992                *output << "ProcPEPsi\t" << config::ProcPEPsi << "\t# for parallel computing: share wave functions" << endl;
     993                *output << "DoOutVis\t" << config::DoOutVis << "\t# Output data for OpenDX" << endl;
     994                *output << "DoOutMes\t" << config::DoOutMes << "\t# Output data for measurements" << endl;
     995                *output << "DoOutOrbitals\t" << config::DoOutOrbitals << "\t# Output all Orbitals" << endl;
     996                *output << "DoOutCurr\t" << config::DoOutCurrent << "\t# Ouput current density for OpenDx" << endl;
     997                *output << "DoOutNICS\t" << config::DoOutNICS << "\t# Output Nucleus independent current shieldings" << endl;
     998                *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
     999                *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
     1000                *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
     1001                *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
     1002                *output << "VectorPlane\t" << config::VectorPlane << "\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot" << endl;
     1003                *output << "VectorCut\t" << config::VectorCut << "\t# Cut plane axis value" << endl;
     1004                *output << "AddGramSch\t" << config::UseAddGramSch << "\t# Additional GramSchmidtOrtogonalization to be safe" << endl;
     1005                *output << "Seed\t\t" << config::Seed << "\t# initial value for random seed for Psi coefficients" << endl;
     1006                *output << endl;
     1007                *output << "MaxOuterStep\t" << config::MaxOuterStep << "\t# number of MolecularDynamics/Structure optimization steps" << endl;
     1008                *output << "Deltat\t" << config::Deltat << "\t# time per MD step" << endl;
     1009                *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
     1010                *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
     1011                *output << "TargetTemp\t" << config::TargetTemp << "\t# Target temperature" << endl;
     1012                *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
     1013                *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
     1014                *output << endl;
     1015                *output << "# Values specifying when to stop" << endl;
     1016                *output << "MaxMinStep\t" << config::MaxMinStep << "\t# Maximum number of steps" << endl;
     1017                *output << "RelEpsTotalE\t" << config::RelEpsTotalEnergy << "\t# relative change in total energy" << endl;
     1018                *output << "RelEpsKineticE\t" << config::RelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
     1019                *output << "MaxMinStopStep\t" << config::MaxMinStopStep << "\t# check every ..th steps" << endl;
     1020                *output << "MaxMinGapStopStep\t" << config::MaxMinGapStopStep << "\t# check every ..th steps" << endl;
     1021                *output << endl;
     1022                *output << "# Values specifying when to stop for INIT, otherwise same as above" << endl;
     1023                *output << "MaxInitMinStep\t" << config::MaxInitMinStep << "\t# Maximum number of steps" << endl;
     1024                *output << "InitRelEpsTotalE\t" << config::InitRelEpsTotalEnergy << "\t# relative change in total energy" << endl;
     1025                *output << "InitRelEpsKineticE\t" << config::InitRelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
     1026                *output << "InitMaxMinStopStep\t" << config::InitMaxMinStopStep << "\t# check every ..th steps" << endl;
     1027                *output << "InitMaxMinGapStopStep\t" << config::InitMaxMinGapStopStep << "\t# check every ..th steps" << endl;
     1028                *output << endl;
     1029                *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
     1030                *output << mol->cell_size[0] << "\t" << endl;
     1031                *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
     1032                *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
     1033                // FIXME
     1034                *output << endl;
     1035                *output << "ECut\t\t" << config::ECut << "\t# energy cutoff for discretization in Hartrees" << endl;
     1036                *output << "MaxLevel\t" << config::MaxLevel << "\t# number of different levels in the code, >=2" << endl;
     1037                *output << "Level0Factor\t" << config::Lev0Factor << "\t# factor by which node number increases from S to 0 level" << endl;
     1038                *output << "RiemannTensor\t" << config::RiemannTensor << "\t# (Use metric)" << endl;
     1039                switch (config::RiemannTensor) {
     1040                        case 0: //UseNoRT
     1041                                break;
     1042                        case 1: // UseRT
     1043                                *output << "RiemannLevel\t" << config::RiemannLevel << "\t# Number of Riemann Levels" << endl;
     1044                                *output << "LevRFactor\t" << config::LevRFactor << "\t# factor by which node number increases from 0 to R level from" << endl;
     1045                                break;
     1046                }
     1047                *output << "PsiType\t\t" << config::PsiType << "\t# 0 - doubly occupied, 1 - SpinUp,SpinDown" << endl;
     1048                // write out both types for easier changing afterwards
     1049        //      switch (PsiType) {
     1050        //              case 0:
     1051                                *output << "MaxPsiDouble\t" << config::MaxPsiDouble << "\t# here: specifying both maximum number of SpinUp- and -Down-states" << endl;
     1052        //                      break;
     1053        //              case 1:
     1054                                *output << "PsiMaxNoUp\t" << config::PsiMaxNoUp << "\t# here: specifying maximum number of SpinUp-states" << endl;
     1055                                *output << "PsiMaxNoDown\t" << config::PsiMaxNoDown << "\t# here: specifying maximum number of SpinDown-states" << endl;
     1056        //                      break;
     1057        //      }
     1058                *output << "AddPsis\t\t" << config::AddPsis << "\t# Additional unoccupied Psis for bandgap determination" << endl;
     1059                *output << endl;
     1060                *output << "RCut\t\t" << config::RCut << "\t# R-cut for the ewald summation" << endl;
     1061                *output << "StructOpt\t" << config::StructOpt << "\t# Do structure optimization beforehand" << endl;
     1062                *output << "IsAngstroem\t" << config::IsAngstroem << "\t# 0 - Bohr, 1 - Angstroem" << endl;
     1063                *output << "RelativeCoord\t" << config::RelativeCoord << "\t# whether ion coordinates are relative (1) or absolute (0)" << endl;
     1064                *output << "MaxTypes\t" << mol->ElementCount << "\t# maximum number of different ion types" << endl;
     1065                *output << endl;
     1066                result = result && mol->Checkout(output);
     1067                if (mol->MDSteps <=1 )
     1068                        result = result && mol->Output(output);
     1069                else
     1070                        result = result && mol->OutputTrajectories(output);
     1071                output->close();
     1072                output->clear();
     1073                delete(output);
     1074                return result;
     1075        } else
     1076                return false;
    10771077};
    10781078
     
    10841084bool config::SaveMPQC(const char *filename, molecule *mol) const
    10851085{
    1086   int ElementNo = 0;
    1087   int AtomNo;
    1088   atom *Walker = NULL;
    1089   element *runner = NULL;
    1090   Vector *center = NULL;
    1091   ofstream *output = NULL;
    1092   stringstream *fname = NULL;
    1093  
    1094   // first without hessian
    1095   fname = new stringstream;
    1096   *fname << filename << ".in";
    1097   output = new ofstream(fname->str().c_str(), ios::out);
    1098   *output << "% Created by MoleCuilder" << endl;
    1099   *output << "mpqc: (" << endl;
    1100   *output << "\tsavestate = no" << endl;
    1101   *output << "\tdo_gradient = yes" << endl;
    1102   *output << "\tmole<CLHF>: (" << endl;
    1103   *output << "\t\tmaxiter = 200" << endl;
    1104   *output << "\t\tbasis = $:basis" << endl;
    1105   *output << "\t\tmolecule = $:molecule" << endl;
    1106   *output << "\t)" << endl;
    1107   *output << ")" << endl;
    1108   *output << "molecule<Molecule>: (" << endl;
    1109   *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
    1110   *output << "\t{ atoms geometry } = {" << endl;
    1111   center = mol->DetermineCenterOfAll(output);
    1112   // output of atoms
    1113   runner = mol->elemente->start;
    1114   while (runner->next != mol->elemente->end) { // go through every element
    1115     runner = runner->next;
    1116     if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
    1117       ElementNo++;
    1118       AtomNo = 0;
    1119       Walker = mol->start;
    1120       while (Walker->next != mol->end) { // go through every atom of this element
    1121         Walker = Walker->next;
    1122         if (Walker->type == runner) { // if this atom fits to element
    1123           AtomNo++;
    1124           *output << "\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
    1125         }
    1126       }
    1127     }
    1128   }
    1129   delete(center);
    1130   *output << "\t}" << endl;
    1131   *output << ")" << endl;
    1132   *output << "basis<GaussianBasisSet>: (" << endl;
    1133   *output << "\tname = \"3-21G\"" << endl;
    1134   *output << "\tmolecule = $:molecule" << endl;
    1135   *output << ")" << endl;
    1136   output->close();
    1137   delete(output);
    1138   delete(fname);
    1139 
    1140   // second with hessian
    1141   fname = new stringstream;
    1142   *fname << filename << ".hess.in";
    1143   output = new ofstream(fname->str().c_str(), ios::out);
    1144   *output << "% Created by MoleCuilder" << endl;
    1145   *output << "mpqc: (" << endl;
    1146   *output << "\tsavestate = no" << endl;
    1147   *output << "\tdo_gradient = yes" << endl;
    1148   *output << "\tmole<CLHF>: (" << endl;
    1149   *output << "\t\tmaxiter = 200" << endl;
    1150   *output << "\t\tbasis = $:basis" << endl;
    1151   *output << "\t\tmolecule = $:molecule" << endl;
    1152   *output << "\t)" << endl;
    1153   *output << "\tfreq<MolecularFrequencies>: (" << endl;
    1154   *output << "\t\tmolecule=$:molecule" << endl;
    1155   *output << "\t)" << endl;
    1156   *output << ")" << endl;
    1157   *output << "molecule<Molecule>: (" << endl;
    1158   *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
    1159   *output << "\t{ atoms geometry } = {" << endl;
    1160   center = mol->DetermineCenterOfAll(output);
    1161   // output of atoms
    1162   runner = mol->elemente->start;
    1163   while (runner->next != mol->elemente->end) { // go through every element
    1164     runner = runner->next;
    1165     if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
    1166       ElementNo++;
    1167       AtomNo = 0;
    1168       Walker = mol->start;
    1169       while (Walker->next != mol->end) { // go through every atom of this element
    1170         Walker = Walker->next;
    1171         if (Walker->type == runner) { // if this atom fits to element
    1172           AtomNo++;
    1173           *output << "\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
    1174         }
    1175       }
    1176     }
    1177   }
    1178   delete(center);
    1179   *output << "\t}" << endl;
    1180   *output << ")" << endl;
    1181   *output << "basis<GaussianBasisSet>: (" << endl;
    1182   *output << "\tname = \"3-21G\"" << endl;
    1183   *output << "\tmolecule = $:molecule" << endl;
    1184   *output << ")" << endl;
    1185   output->close();
    1186   delete(output);
    1187   delete(fname);
    1188  
    1189   return true;
     1086        int ElementNo = 0;
     1087        int AtomNo;
     1088        atom *Walker = NULL;
     1089        element *runner = NULL;
     1090        Vector *center = NULL;
     1091        ofstream *output = NULL;
     1092        stringstream *fname = NULL;
     1093       
     1094        // first without hessian
     1095        fname = new stringstream;
     1096        *fname << filename << ".in";
     1097        output = new ofstream(fname->str().c_str(), ios::out);
     1098        *output << "% Created by MoleCuilder" << endl;
     1099        *output << "mpqc: (" << endl;
     1100        *output << "\tsavestate = no" << endl;
     1101        *output << "\tdo_gradient = yes" << endl;
     1102        *output << "\tmole<MBPT2>: (" << endl;
     1103        *output << "\t\tmaxiter = 200" << endl;
     1104        *output << "\t\tbasis = $:basis" << endl;
     1105        *output << "\t\tmolecule = $:molecule" << endl;
     1106        *output << "\t\treference<CLHF>: (" << endl;
     1107        *output << "\t\t\tbasis = $:basis" << endl;
     1108        *output << "\t\t\tmolecule = $:molecule" << endl;
     1109        *output << "\t\t)" << endl;
     1110        *output << "\t)" << endl;
     1111        *output << ")" << endl;
     1112        *output << "molecule<Molecule>: (" << endl;
     1113        *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
     1114        *output << "\t{ atoms geometry } = {" << endl;
     1115        center = mol->DetermineCenterOfAll(output);
     1116        // output of atoms
     1117        runner = mol->elemente->start;
     1118        while (runner->next != mol->elemente->end) { // go through every element
     1119                runner = runner->next;
     1120                if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
     1121                        ElementNo++;
     1122                        AtomNo = 0;
     1123                        Walker = mol->start;
     1124                        while (Walker->next != mol->end) { // go through every atom of this element
     1125                                Walker = Walker->next;
     1126                                if (Walker->type == runner) { // if this atom fits to element
     1127                                        AtomNo++;
     1128                                        *output << "\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
     1129                                }
     1130                        }
     1131                }
     1132        }
     1133        delete(center);
     1134        *output << "\t}" << endl;
     1135        *output << ")" << endl;
     1136        *output << "basis<GaussianBasisSet>: (" << endl;
     1137        *output << "\tname = \"3-21G\"" << endl;
     1138        *output << "\tmolecule = $:molecule" << endl;
     1139        *output << ")" << endl;
     1140        output->close();
     1141        delete(output);
     1142        delete(fname);
     1143
     1144        // second with hessian
     1145        fname = new stringstream;
     1146        *fname << filename << ".hess.in";
     1147        output = new ofstream(fname->str().c_str(), ios::out);
     1148        *output << "% Created by MoleCuilder" << endl;
     1149        *output << "mpqc: (" << endl;
     1150        *output << "\tsavestate = no" << endl;
     1151        *output << "\tdo_gradient = yes" << endl;
     1152        *output << "\tmole<CLHF>: (" << endl;
     1153        *output << "\t\tmaxiter = 200" << endl;
     1154        *output << "\t\tbasis = $:basis" << endl;
     1155        *output << "\t\tmolecule = $:molecule" << endl;
     1156        *output << "\t)" << endl;
     1157        *output << "\tfreq<MolecularFrequencies>: (" << endl;
     1158        *output << "\t\tmolecule=$:molecule" << endl;
     1159        *output << "\t)" << endl;
     1160        *output << ")" << endl;
     1161        *output << "molecule<Molecule>: (" << endl;
     1162        *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
     1163        *output << "\t{ atoms geometry } = {" << endl;
     1164        center = mol->DetermineCenterOfAll(output);
     1165        // output of atoms
     1166        runner = mol->elemente->start;
     1167        while (runner->next != mol->elemente->end) { // go through every element
     1168                runner = runner->next;
     1169                if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
     1170                        ElementNo++;
     1171                        AtomNo = 0;
     1172                        Walker = mol->start;
     1173                        while (Walker->next != mol->end) { // go through every atom of this element
     1174                                Walker = Walker->next;
     1175                                if (Walker->type == runner) { // if this atom fits to element
     1176                                        AtomNo++;
     1177                                        *output << "\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
     1178                                }
     1179                        }
     1180                }
     1181        }
     1182        delete(center);
     1183        *output << "\t}" << endl;
     1184        *output << ")" << endl;
     1185        *output << "basis<GaussianBasisSet>: (" << endl;
     1186        *output << "\tname = \"3-21G\"" << endl;
     1187        *output << "\tmolecule = $:molecule" << endl;
     1188        *output << ")" << endl;
     1189        output->close();
     1190        delete(output);
     1191        delete(fname);
     1192       
     1193        return true;
    11901194};
    11911195
     
    11991203 * \param name Name of value in file (at least 3 chars!)
    12001204 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
    1201  *        (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
    1202  *        best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
    1203  *        counted from this unresetted position!)
     1205 *                              (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
     1206 *                              best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
     1207 *                              counted from this unresetted position!)
    12041208 * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
    12051209 * \param yth In grid case specifying column number, otherwise the yth \a name matching line
     
    12121216 */
    12131217int config::ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical) {
    1214   int i,j;  // loop variables
    1215   int length = 0, maxlength = -1;
    1216   long file_position = file->tellg(); // mark current position
    1217   char *dummy1, *dummy, *free_dummy;    // pointers in the line that is read in per step
    1218   dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char), "config::ParseForParameter: *free_dummy");
    1219 
    1220   //fprintf(stderr,"Parsing for %s\n",name); 
    1221   if (repetition == 0)
    1222     //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
    1223     return 0;
    1224 
    1225   int line = 0; // marks line where parameter was found
    1226   int found = (type >= grid) ? 0 : (-yth + 1);  // marks if yth parameter name was found
    1227   while((found != repetition)) {
    1228     dummy1 = dummy = free_dummy;
    1229     do {
    1230       file->getline(dummy1, 256); // Read the whole line
    1231       if (file->eof()) {
    1232         if ((critical) && (found == 0)) {
    1233           Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1234           //Error(InitReading, name);
    1235           fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
    1236           exit(255);
    1237         } else {
    1238           //if (!sequential)
    1239           file->clear();
    1240           file->seekg(file_position, ios::beg);  // rewind to start position
    1241           Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");         
    1242           return 0;
    1243         }
    1244       }
    1245       line++;
    1246     } while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
    1247    
    1248     // C++ getline removes newline at end, thus re-add
    1249     if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
    1250       i = strlen(dummy1);
    1251       dummy1[i] = '\n';
    1252       dummy1[i+1] = '\0';
    1253     }
    1254     //fprintf(stderr,"line %i ends at %i, newline at %i\n", line, strlen(dummy1), strchr(dummy1,'\n')-free_dummy);
    1255 
    1256     if (dummy1 == NULL) {
    1257       if (verbose) fprintf(stderr,"Error reading line %i\n",line);
    1258     } else {
    1259       //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
    1260     }
    1261     // Seek for possible end of keyword on line if given ...
    1262     if (name != NULL) {
    1263       dummy = strchr(dummy1,'\t');  // set dummy on first tab or space which ever's nearer
    1264       if (dummy == NULL) {
    1265         dummy = strchr(dummy1, ' ');  // if not found seek for space
    1266         while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
    1267           dummy++;
    1268       }
    1269       if (dummy == NULL) {
    1270         dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
    1271         //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
    1272         //Free((void **)&free_dummy);
    1273         //Error(FileOpenParams, NULL);     
    1274       } else {
    1275         //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
    1276       }
    1277     } else dummy = dummy1;
    1278     // ... and check if it is the keyword!
    1279     //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
    1280     if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
    1281       found++; // found the parameter!
    1282       //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);   
    1283      
    1284       if (found == repetition) {
    1285         for (i=0;i<xth;i++) { // i = rows
    1286           if (type >= grid) {
    1287             // grid structure means that grid starts on the next line, not right after keyword
    1288             dummy1 = dummy = free_dummy;
    1289             do {
    1290               file->getline(dummy1, 256); // Read the whole line, skip commentary and empty ones
    1291               if (file->eof()) {
    1292                 if ((critical) && (found == 0)) {
    1293                   Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1294                   //Error(InitReading, name);
    1295                   fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
    1296                   exit(255);
    1297                 } else {
    1298                   //if (!sequential)
    1299                   file->clear();
    1300                   file->seekg(file_position, ios::beg);  // rewind to start position
    1301                   Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1302                   return 0;
    1303                 }
    1304               }
    1305               line++;
    1306             } while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
    1307             if (dummy1 == NULL){
    1308               if (verbose) fprintf(stderr,"Error reading line %i\n", line);
    1309             } else {
    1310               //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
    1311             }
    1312           } else { // simple int, strings or doubles start in the same line
    1313             while ((*dummy == '\t') || (*dummy == ' '))  // skip interjacent tabs and spaces
    1314               dummy++;
    1315           }
    1316           // C++ getline removes newline at end, thus re-add
    1317           if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
    1318             j = strlen(dummy1);
    1319             dummy1[j] = '\n';
    1320             dummy1[j+1] = '\0';
    1321           }
    1322  
    1323           int start = (type >= grid) ? 0 : yth-1 ;
    1324           for (j=start;j<yth;j++) { // j = columns
    1325             // check for lower triangular area and upper triangular area
    1326             if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
    1327               *((double *)value) = 0.0;
    1328               fprintf(stderr,"%f\t",*((double *)value));
    1329               value = (void *)((long)value + sizeof(double));
    1330               //value += sizeof(double);
    1331             } else {
    1332               // otherwise we must skip all interjacent tabs and spaces and find next value
    1333               dummy1 = dummy;
    1334               dummy = strchr(dummy1, '\t'); // seek for tab or space
    1335               if (dummy == NULL)
    1336                 dummy = strchr(dummy1, ' ');  // if not found seek for space
    1337               if (dummy == NULL) { // if still zero returned ...
    1338                 dummy = strchr(dummy1, '\n');  // ... at line end then
    1339                 if ((j < yth-1) && (type < 4)) {  // check if xth value or not yet
    1340                   if (critical) {
    1341                     if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
    1342                     Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1343                     //return 0;
    1344                     exit(255);
    1345                     //Error(FileOpenParams, NULL);     
    1346                   } else {
    1347                     //if (!sequential)
    1348                     file->clear();
    1349                     file->seekg(file_position, ios::beg);  // rewind to start position
    1350                     Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1351                     return 0;
    1352                   }
    1353                 }
    1354               } else {
    1355                 //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
    1356               }
    1357               if (*dummy1 == '#') {
    1358                 // found comment, skipping rest of line
    1359                 //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
    1360                 if (!sequential) { // here we need it!
    1361                   file->seekg(file_position, ios::beg);  // rewind to start position
    1362                 }
    1363                 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1364                 return 0;
    1365               }
    1366               //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
    1367               switch(type) {
    1368                 case (row_int):
    1369                   *((int *)value) = atoi(dummy1);
    1370                   if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
    1371                   if (verbose) fprintf(stderr,"%i\t",*((int *)value));
    1372                     value = (void *)((long)value + sizeof(int));
    1373                     //value += sizeof(int);
    1374                   break;
    1375                 case(row_double):
    1376                 case(grid):
    1377                 case(lower_trigrid):
    1378                 case(upper_trigrid):
    1379                   *((double *)value) = atof(dummy1);
    1380                   if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
    1381                   if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
    1382                   value = (void *)((long)value + sizeof(double));
    1383                   //value += sizeof(double);
    1384                   break;
    1385                 case(double_type):
    1386                   *((double *)value) = atof(dummy1);
    1387                   if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
    1388                   //value += sizeof(double);
    1389                   break;
    1390                 case(int_type):
    1391                   *((int *)value) = atoi(dummy1);
    1392                   if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
    1393                   //value += sizeof(int);
    1394                   break;
    1395                 default:
    1396                 case(string_type):
    1397                   if (value != NULL) {
    1398                     //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
    1399                     maxlength = MAXSTRINGSIZE;
    1400                     length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
    1401                     strncpy((char *)value, dummy1, length);  // copy as much
    1402                     ((char *)value)[length] = '\0';  // and set end marker
    1403                     if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
    1404                     //value += sizeof(char);
    1405                   } else {
    1406                   }
    1407                 break;
    1408               }
    1409             }
    1410             while (*dummy == '\t')
    1411               dummy++;
    1412           }
    1413         }
    1414       }
    1415     }
    1416   }
    1417   if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
    1418   Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    1419   if (!sequential) {
    1420     file->clear();
    1421     file->seekg(file_position, ios::beg);  // rewind to start position
    1422   }
    1423   //fprintf(stderr, "End of Parsing\n\n");
    1424  
    1425   return (found); // true if found, false if not
     1218        int i,j;        // loop variables
     1219        int length = 0, maxlength = -1;
     1220        long file_position = file->tellg(); // mark current position
     1221        char *dummy1, *dummy, *free_dummy;              // pointers in the line that is read in per step
     1222        dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char), "config::ParseForParameter: *free_dummy");
     1223
     1224        //fprintf(stderr,"Parsing for %s\n",name);     
     1225        if (repetition == 0)
     1226                //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
     1227                return 0;
     1228
     1229        int line = 0; // marks line where parameter was found
     1230        int found = (type >= grid) ? 0 : (-yth + 1);    // marks if yth parameter name was found
     1231        while((found != repetition)) {
     1232                dummy1 = dummy = free_dummy;
     1233                do {
     1234                        file->getline(dummy1, 256); // Read the whole line
     1235                        if (file->eof()) {
     1236                                if ((critical) && (found == 0)) {
     1237                                        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1238                                        //Error(InitReading, name);
     1239                                        fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
     1240                                        exit(255);
     1241                                } else {
     1242                                        //if (!sequential)
     1243                                        file->clear();
     1244                                        file->seekg(file_position, ios::beg);   // rewind to start position
     1245                                        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");                                   
     1246                                        return 0;
     1247                                }
     1248                        }
     1249                        line++;
     1250                } while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
     1251               
     1252                // C++ getline removes newline at end, thus re-add
     1253                if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
     1254                        i = strlen(dummy1);
     1255                        dummy1[i] = '\n';
     1256                        dummy1[i+1] = '\0';
     1257                }
     1258                //fprintf(stderr,"line %i ends at %i, newline at %i\n", line, strlen(dummy1), strchr(dummy1,'\n')-free_dummy);
     1259
     1260                if (dummy1 == NULL) {
     1261                        if (verbose) fprintf(stderr,"Error reading line %i\n",line);
     1262                } else {
     1263                        //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
     1264                }
     1265                // Seek for possible end of keyword on line if given ...
     1266                if (name != NULL) {
     1267                        dummy = strchr(dummy1,'\t');    // set dummy on first tab or space which ever's nearer
     1268                        if (dummy == NULL) {
     1269                                dummy = strchr(dummy1, ' ');    // if not found seek for space
     1270                                while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))                // skip some more tabs and spaces if necessary
     1271                                        dummy++;
     1272                        }
     1273                        if (dummy == NULL) {
     1274                                dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
     1275                                //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
     1276                                //Free((void **)&free_dummy);
     1277                                //Error(FileOpenParams, NULL);                 
     1278                        } else {
     1279                                //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
     1280                        }
     1281                } else dummy = dummy1;
     1282                // ... and check if it is the keyword!
     1283                //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
     1284                if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
     1285                        found++; // found the parameter!
     1286                        //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);         
     1287                       
     1288                        if (found == repetition) {
     1289                                for (i=0;i<xth;i++) { // i = rows
     1290                                        if (type >= grid) {
     1291                                                // grid structure means that grid starts on the next line, not right after keyword
     1292                                                dummy1 = dummy = free_dummy;
     1293                                                do {
     1294                                                        file->getline(dummy1, 256); // Read the whole line, skip commentary and empty ones
     1295                                                        if (file->eof()) {
     1296                                                                if ((critical) && (found == 0)) {
     1297                                                                        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1298                                                                        //Error(InitReading, name);
     1299                                                                        fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
     1300                                                                        exit(255);
     1301                                                                } else {
     1302                                                                        //if (!sequential)
     1303                                                                        file->clear();
     1304                                                                        file->seekg(file_position, ios::beg);   // rewind to start position
     1305                                                                        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1306                                                                        return 0;
     1307                                                                }
     1308                                                        }
     1309                                                        line++;
     1310                                                } while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
     1311                                                if (dummy1 == NULL){
     1312                                                        if (verbose) fprintf(stderr,"Error reading line %i\n", line);
     1313                                                } else {
     1314                                                        //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
     1315                                                }
     1316                                        } else { // simple int, strings or doubles start in the same line
     1317                                                while ((*dummy == '\t') || (*dummy == ' '))      // skip interjacent tabs and spaces
     1318                                                        dummy++;
     1319                                        }
     1320                                        // C++ getline removes newline at end, thus re-add
     1321                                        if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
     1322                                                j = strlen(dummy1);
     1323                                                dummy1[j] = '\n';
     1324                                                dummy1[j+1] = '\0';
     1325                                        }
     1326       
     1327                                        int start = (type >= grid) ? 0 : yth-1 ;
     1328                                        for (j=start;j<yth;j++) { // j = columns
     1329                                                // check for lower triangular area and upper triangular area
     1330                                                if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
     1331                                                        *((double *)value) = 0.0;
     1332                                                        fprintf(stderr,"%f\t",*((double *)value));
     1333                                                        value = (void *)((long)value + sizeof(double));
     1334                                                        //value += sizeof(double);
     1335                                                } else {
     1336                                                        // otherwise we must skip all interjacent tabs and spaces and find next value
     1337                                                        dummy1 = dummy;
     1338                                                        dummy = strchr(dummy1, '\t'); // seek for tab or space
     1339                                                        if (dummy == NULL)
     1340                                                                dummy = strchr(dummy1, ' ');    // if not found seek for space
     1341                                                        if (dummy == NULL) { // if still zero returned ...
     1342                                                                dummy = strchr(dummy1, '\n');   // ... at line end then
     1343                                                                if ((j < yth-1) && (type < 4)) {        // check if xth value or not yet
     1344                                                                        if (critical) {
     1345                                                                                if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
     1346                                                                                Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1347                                                                                //return 0;
     1348                                                                                exit(255);
     1349                                                                                //Error(FileOpenParams, NULL);                 
     1350                                                                        } else {
     1351                                                                                //if (!sequential)
     1352                                                                                file->clear();
     1353                                                                                file->seekg(file_position, ios::beg);   // rewind to start position
     1354                                                                                Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1355                                                                                return 0;
     1356                                                                        }
     1357                                                                }
     1358                                                        } else {
     1359                                                                //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
     1360                                                        }
     1361                                                        if (*dummy1 == '#') {
     1362                                                                // found comment, skipping rest of line
     1363                                                                //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
     1364                                                                if (!sequential) { // here we need it!
     1365                                                                        file->seekg(file_position, ios::beg);   // rewind to start position
     1366                                                                }
     1367                                                                Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1368                                                                return 0;
     1369                                                        }
     1370                                                        //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
     1371                                                        switch(type) {
     1372                                                                case (row_int):
     1373                                                                        *((int *)value) = atoi(dummy1);
     1374                                                                        if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
     1375                                                                        if (verbose) fprintf(stderr,"%i\t",*((int *)value));
     1376                                                                                value = (void *)((long)value + sizeof(int));
     1377                                                                                //value += sizeof(int);
     1378                                                                        break;
     1379                                                                case(row_double):
     1380                                                                case(grid):
     1381                                                                case(lower_trigrid):
     1382                                                                case(upper_trigrid):
     1383                                                                        *((double *)value) = atof(dummy1);
     1384                                                                        if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
     1385                                                                        if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
     1386                                                                        value = (void *)((long)value + sizeof(double));
     1387                                                                        //value += sizeof(double);
     1388                                                                        break;
     1389                                                                case(double_type):
     1390                                                                        *((double *)value) = atof(dummy1);
     1391                                                                        if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
     1392                                                                        //value += sizeof(double);
     1393                                                                        break;
     1394                                                                case(int_type):
     1395                                                                        *((int *)value) = atoi(dummy1);
     1396                                                                        if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
     1397                                                                        //value += sizeof(int);
     1398                                                                        break;
     1399                                                                default:
     1400                                                                case(string_type):
     1401                                                                        if (value != NULL) {
     1402                                                                                //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
     1403                                                                                maxlength = MAXSTRINGSIZE;
     1404                                                                                length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
     1405                                                                                strncpy((char *)value, dummy1, length); // copy as much
     1406                                                                                ((char *)value)[length] = '\0'; // and set end marker
     1407                                                                                if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
     1408                                                                                //value += sizeof(char);
     1409                                                                        } else {
     1410                                                                        }
     1411                                                                break;
     1412                                                        }
     1413                                                }
     1414                                                while (*dummy == '\t')
     1415                                                        dummy++;
     1416                                        }
     1417                                }
     1418                        }
     1419                }
     1420        }
     1421        if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
     1422        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1423        if (!sequential) {
     1424                file->clear();
     1425                file->seekg(file_position, ios::beg);   // rewind to start position
     1426        }
     1427        //fprintf(stderr, "End of Parsing\n\n");
     1428       
     1429        return (found); // true if found, false if not
    14261430}
Note: See TracChangeset for help on using the changeset viewer.