source: src/config.cpp@ fa649a

/** \file config.cpp
 *
 * Function implementations for the class config.
 *
 */

#include "atom.hpp"
#include "config.hpp"
#include "element.hpp"
#include "helpers.hpp"
#include "lists.hpp"
#include "molecule.hpp"
#include "memoryallocator.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"

/******************************** Functions for class ConfigFileBuffer **********************/

/** Structure containing compare function for Ion_Type sorting.
 */
struct IonTypeCompare {
  bool operator()(const char* s1, const char *s2) const {
    char number1[8];
    char number2[8];
    char *dummy1, *dummy2;
    //cout << s1 << "  " << s2 << endl;
    dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
    dummy2 = strchr(dummy1, '_');
    strncpy(number1, dummy1, dummy2-dummy1); // copy the number
    number1[dummy2-dummy1]='\0';
    dummy1 = strchr(s2, '_')+sizeof(char)*5;  // go just after "Ion_Type"
    dummy2 = strchr(dummy1, '_');
    strncpy(number2, dummy1, dummy2-dummy1); // copy the number
    number2[dummy2-dummy1]='\0';
    if (atoi(number1) != atoi(number2))
      return (atoi(number1) < atoi(number2));
    else {
      dummy1 = strchr(s1, '_')+sizeof(char);
      dummy1 = strchr(dummy1, '_')+sizeof(char);
      dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
      strncpy(number1, dummy1, dummy2-dummy1); // copy the number
      number1[dummy2-dummy1]='\0';
      dummy1 = strchr(s2, '_')+sizeof(char);
      dummy1 = strchr(dummy1, '_')+sizeof(char);
      dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
      strncpy(number2, dummy1, dummy2-dummy1); // copy the number
      number2[dummy2-dummy1]='\0';
      return (atoi(number1) < atoi(number2));
    }
  }
};

/** Constructor for ConfigFileBuffer class.
 */
ConfigFileBuffer::ConfigFileBuffer() : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
{
};

/** Constructor for ConfigFileBuffer class with filename to be parsed.
 * \param *filename file name
 */
ConfigFileBuffer::ConfigFileBuffer(const char * const filename) : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
{
  ifstream *file = NULL;
  char line[MAXSTRINGSIZE];

  // prescan number of lines
  file= new ifstream(filename);
  if (file == NULL) {
    cerr << "ERROR: config file " << filename << " missing!" << endl;
    return;
  }
  NoLines = 0; // we're overcounting by one
  long file_position = file->tellg(); // mark current position
  do {
    file->getline(line, 256);
    NoLines++;
  } while (!file->eof());
  file->clear();
  file->seekg(file_position, ios::beg);
  cout << Verbose(1) << NoLines-1 << " lines were recognized." << endl;

  // allocate buffer's 1st dimension
  if (buffer != NULL) {
    cerr << "ERROR: FileBuffer->buffer is not NULL!" << endl;
    return;
  } else
    buffer = Malloc<char*>(NoLines, "ConfigFileBuffer::ConfigFileBuffer: **buffer");

  // scan each line and put into buffer
  int lines=0;
  int i;
  do {
    buffer[lines] = Malloc<char>(MAXSTRINGSIZE, "ConfigFileBuffer::ConfigFileBuffer: *buffer[]");
    file->getline(buffer[lines], MAXSTRINGSIZE-1);
    i = strlen(buffer[lines]);
    buffer[lines][i] = '\n';
    buffer[lines][i+1] = '\0';
    lines++;
  } while((!file->eof()) && (lines < NoLines));
  cout << Verbose(1) << lines-1 << " lines were read into the buffer." << endl;

  // close and exit
  file->close();
  file->clear();
  delete(file);
}

/** Destructor for ConfigFileBuffer class.
 */
ConfigFileBuffer::~ConfigFileBuffer()
{
  for(int i=0;i<NoLines;++i)
    Free(&buffer[i]);
  Free(&buffer);
  Free(&LineMapping);
}


/** Create trivial mapping.
 */
void ConfigFileBuffer::InitMapping()
{
  LineMapping = Malloc<int>(NoLines, "ConfigFileBuffer::InitMapping: *LineMapping");
  for (int i=0;i<NoLines;i++)
    LineMapping[i] = i;
}

/** Creates a mapping for the \a *FileBuffer's lines containing the Ion_Type keyword such that they are sorted.
 * \a *map on return contains a list of NoAtom entries such that going through the list, yields indices to the
 * lines in \a *FileBuffer in a sorted manner of the Ion_Type?_? keywords. We assume that ConfigFileBuffer::CurrentLine
 * points to first Ion_Type entry.
 * \param *FileBuffer pointer to buffer structure
 * \param NoAtoms of subsequent lines to look at
 */
void ConfigFileBuffer::MapIonTypesInBuffer(const int NoAtoms)
{
  map<const char *, int, IonTypeCompare> LineList;
  if (LineMapping == NULL) {
    cerr << "map pointer is NULL: " << LineMapping << endl;
    return;
  }

  // put all into hashed map
  for (int i=0; i<NoAtoms; ++i) {
    LineList.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i));
  }

  // fill map
  int nr=0;
  for (map<const char *, int, IonTypeCompare>::iterator runner = LineList.begin(); runner != LineList.end(); ++runner) {
    if (CurrentLine+nr < NoLines)
      LineMapping[CurrentLine+(nr++)] = runner->second;
    else
      cerr << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl;
  }
}

/************************************* Functions for class config ***************************/

/** Constructor for config file class.
 */
config::config() : BG(NULL), PsiType(0), MaxPsiDouble(0), PsiMaxNoUp(0), PsiMaxNoDown(0), MaxMinStopStep(1), InitMaxMinStopStep(1), ProcPEGamma(8), ProcPEPsi(1), configpath(NULL),
    configname(NULL), FastParsing(false), Deltat(0.01), basis(""), databasepath(NULL), DoConstrainedMD(0), MaxOuterStep(0), Thermostat(4), ThermostatImplemented(NULL),
    ThermostatNames(NULL), TempFrequency(2.5), alpha(0.), HooverMass(0.), TargetTemp(0.00095004455), ScaleTempStep(25),  mainname(NULL), defaultpath(NULL), pseudopotpath(NULL),
    DoOutVis(0), DoOutMes(1), DoOutNICS(0), DoOutOrbitals(0), DoOutCurrent(0), DoFullCurrent(0), DoPerturbation(0), DoWannier(0), CommonWannier(0), SawtoothStart(0.01),
    VectorPlane(0), VectorCut(0.), UseAddGramSch(1), Seed(1), OutVisStep(10), OutSrcStep(5), MaxPsiStep(0), EpsWannier(1e-7), MaxMinStep(100), RelEpsTotalEnergy(1e-7),
    RelEpsKineticEnergy(1e-5), MaxMinGapStopStep(0), MaxInitMinStep(100), InitRelEpsTotalEnergy(1e-5), InitRelEpsKineticEnergy(1e-4), InitMaxMinGapStopStep(0), ECut(128.),
    MaxLevel(5), RiemannTensor(0), LevRFactor(0), RiemannLevel(0), Lev0Factor(2), RTActualUse(0), AddPsis(0), RCut(20.), StructOpt(0), IsAngstroem(1), RelativeCoord(0),
    MaxTypes(0) {
  mainname = Malloc<char>(MAXSTRINGSIZE,"config constructor: mainname");
  defaultpath = Malloc<char>(MAXSTRINGSIZE,"config constructor: defaultpath");
  pseudopotpath = Malloc<char>(MAXSTRINGSIZE,"config constructor: pseudopotpath");
  databasepath = Malloc<char>(MAXSTRINGSIZE,"config constructor: databasepath");
  configpath = Malloc<char>(MAXSTRINGSIZE,"config constructor: configpath");
  configname = Malloc<char>(MAXSTRINGSIZE,"config constructor: configname");
  strcpy(mainname,"pcp");
  strcpy(defaultpath,"not specified");
  strcpy(pseudopotpath,"not specified");
  configpath[0]='\0';
  configname[0]='\0';
  basis = "3-21G";

  InitThermostats();
};

/** Destructor for config file class.
 */
config::~config()
{
  Free(&mainname);
  Free(&defaultpath);
  Free(&pseudopotpath);
  Free(&databasepath);
  Free(&configpath);
  Free(&configname);
  Free(&ThermostatImplemented);
  for (int j=0;j<MaxThermostats;j++)
    Free(&ThermostatNames[j]);
  Free(&ThermostatNames);
};

/** Initialises variables in class config for Thermostats.
 */
void config::InitThermostats()
{
  ThermostatImplemented = Malloc<int>(MaxThermostats, "config constructor: *ThermostatImplemented");
  ThermostatNames = Malloc<char*>(MaxThermostats, "config constructor: *ThermostatNames");
  for (int j=0;j<MaxThermostats;j++)
    ThermostatNames[j] = Malloc<char>(12, "config constructor: ThermostatNames[]");

  strcpy(ThermostatNames[0],"None");
  ThermostatImplemented[0] = 1;
  strcpy(ThermostatNames[1],"Woodcock");
  ThermostatImplemented[1] = 1;
  strcpy(ThermostatNames[2],"Gaussian");
  ThermostatImplemented[2] = 1;
  strcpy(ThermostatNames[3],"Langevin");
  ThermostatImplemented[3] = 1;
  strcpy(ThermostatNames[4],"Berendsen");
  ThermostatImplemented[4] = 1;
  strcpy(ThermostatNames[5],"NoseHoover");
  ThermostatImplemented[5] = 1;
};

/** Readin of Thermostat related values from parameter file.
 * \param *fb file buffer containing the config file
 */
void config::ParseThermostats(class ConfigFileBuffer * const fb)
{
  char * const thermo = Malloc<char>(12, "IonsInitRead: thermo");
  const int verbose = 0;

  // read desired Thermostat from file along with needed additional parameters
  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
      if (ThermostatImplemented[0] == 1) {
        Thermostat = None;
      } else {
        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
        Thermostat = None;
      }
    } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
      if (ThermostatImplemented[1] == 1) {
        Thermostat = Woodcock;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
      } else {
        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
        Thermostat = None;
      }
    } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
      if (ThermostatImplemented[2] == 1) {
        Thermostat = Gaussian;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
      } else {
        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
        Thermostat = None;
      }
    } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
      if (ThermostatImplemented[3] == 1) {
        Thermostat = Langevin;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
          cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
        } else {
          alpha = 1.;
        }
      } else {
        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
        Thermostat = None;
      }
    } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
      if (ThermostatImplemented[4] == 1) {
        Thermostat = Berendsen;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
      } else {
        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
        Thermostat = None;
      }
    } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
      if (ThermostatImplemented[5] == 1) {
        Thermostat = NoseHoover;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
        alpha = 0.;
      } else {
        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
        Thermostat = None;
      }
    } else {
      cout << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
      Thermostat = None;
    }
  } else {
    if ((MaxOuterStep > 0) && (TargetTemp != 0)) 
      cout << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
    Thermostat = None;
  }
  Free(thermo);
};


/** Displays menu for editing each entry of the config file.
 * Nothing fancy here, just lots of cout << Verbose(0)s for the menu and a switch/case
 * for each entry of the config file structure.
 */
void config::Edit()
{
  char choice;

  do {
    cout << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
    cout << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl;
    cout << Verbose(0) << " B - Default path (for runtime files)" << endl;
    cout << Verbose(0) << " C - Path of pseudopotential files" << endl;
    cout << Verbose(0) << " D - Number of coefficient sharing processes" << endl;
    cout << Verbose(0) << " E - Number of wave function sharing processes" << endl;
    cout << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl;
    cout << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl;
    cout << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl;
    cout << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl;
    cout << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl;
    cout << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl;
    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;
    cout << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl;
    cout << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl;
    cout << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl;
    cout << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl;
    cout << Verbose(0) << " Q - Initial integer value of random number generator" << endl;
    cout << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl;
    cout << Verbose(0) << " T - Output visual after ...th step" << endl;
    cout << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl;
    cout << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl;
    cout << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl;
    cout << Verbose(0) << " Z - Maximum number of minimization iterations" << endl;
    cout << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl;
    cout << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl;
    cout << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl;
    cout << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl;
    cout << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl;
    cout << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
    cout << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
//    cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
    cout << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
    cout << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
    cout << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl;
    cout << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl;
    cout << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl;
    cout << Verbose(0) << " p - Number of Riemann levels" << endl;
    cout << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl;
    cout << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl;
    cout << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl;
    cout << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl;
    cout << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl;
    cout << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl;
    cout << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl;
    cout << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl;
    cout << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl;
    cout << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl;
    cout << Verbose(0) << "=========================================================" << endl;
    cout << Verbose(0) << "INPUT: ";
    cin >> choice;

    switch (choice) {
        case 'A': // mainname
          cout << Verbose(0) << "Old: " << config::mainname << "\t new: ";
          cin >> config::mainname;
          break;
        case 'B': // defaultpath
          cout << Verbose(0) << "Old: " << config::defaultpath << "\t new: ";
          cin >> config::defaultpath;
          break;
        case 'C': // pseudopotpath
          cout << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ";
          cin >> config::pseudopotpath;
          break;

        case 'D': // ProcPEGamma
          cout << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ";
          cin >> config::ProcPEGamma;
          break;
        case 'E': // ProcPEPsi
          cout << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ";
          cin >> config::ProcPEPsi;
          break;
        case 'F': // DoOutVis
          cout << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ";
          cin >> config::DoOutVis;
          break;
        case 'G': // DoOutMes
          cout << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ";
          cin >> config::DoOutMes;
          break;
        case 'H': // DoOutOrbitals
          cout << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ";
          cin >> config::DoOutOrbitals;
          break;
        case 'I': // DoOutCurrent
          cout << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ";
          cin >> config::DoOutCurrent;
          break;
        case 'J': // DoFullCurrent
          cout << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ";
          cin >> config::DoFullCurrent;
          break;
        case 'K': // DoPerturbation
          cout << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ";
          cin >> config::DoPerturbation;
          break;
        case 'L': // CommonWannier
          cout << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ";
          cin >> config::CommonWannier;
          break;
        case 'M': // SawtoothStart
          cout << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ";
          cin >> config::SawtoothStart;
          break;
        case 'N': // VectorPlane
          cout << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ";
          cin >> config::VectorPlane;
          break;
        case 'O': // VectorCut
          cout << Verbose(0) << "Old: " << config::VectorCut << "\t new: ";
          cin >> config::VectorCut;
          break;
        case 'P': // UseAddGramSch
          cout << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ";
          cin >> config::UseAddGramSch;
          break;
        case 'Q': // Seed
          cout << Verbose(0) << "Old: " << config::Seed << "\t new: ";
          cin >> config::Seed;
          break;

        case 'R': // MaxOuterStep
          cout << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ";
          cin >> config::MaxOuterStep;
          break;
        case 'T': // OutVisStep
          cout << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ";
          cin >> config::OutVisStep;
          break;
        case 'U': // OutSrcStep
          cout << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ";
          cin >> config::OutSrcStep;
          break;
        case 'X': // MaxPsiStep
          cout << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ";
          cin >> config::MaxPsiStep;
          break;
        case 'Y': // EpsWannier
          cout << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ";
          cin >> config::EpsWannier;
          break;

        case 'Z': // MaxMinStep
          cout << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ";
          cin >> config::MaxMinStep;
          break;
        case 'a': // RelEpsTotalEnergy
          cout << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ";
          cin >> config::RelEpsTotalEnergy;
          break;
        case 'b': // RelEpsKineticEnergy
          cout << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ";
          cin >> config::RelEpsKineticEnergy;
          break;
        case 'c': // MaxMinStopStep
          cout << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ";
          cin >> config::MaxMinStopStep;
          break;
        case 'e': // MaxInitMinStep
          cout << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ";
          cin >> config::MaxInitMinStep;
          break;
        case 'f': // InitRelEpsTotalEnergy
          cout << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ";
          cin >> config::InitRelEpsTotalEnergy;
          break;
        case 'g': // InitRelEpsKineticEnergy
          cout << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ";
          cin >> config::InitRelEpsKineticEnergy;
          break;
        case 'h': // InitMaxMinStopStep
          cout << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ";
          cin >> config::InitMaxMinStopStep;
          break;

//        case 'j': // BoxLength
//          cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
//          for (int i=0;i<6;i++) {
//            cout << Verbose(0) << "Cell size" << i << ": ";
//            cin >> mol->cell_size[i];
//          }
//          break;

        case 'k': // ECut
          cout << Verbose(0) << "Old: " << config::ECut << "\t new: ";
          cin >> config::ECut;
          break;
        case 'l': // MaxLevel
          cout << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ";
          cin >> config::MaxLevel;
          break;
        case 'm': // RiemannTensor
          cout << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ";
          cin >> config::RiemannTensor;
          break;
        case 'n': // LevRFactor
          cout << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ";
          cin >> config::LevRFactor;
          break;
        case 'o': // RiemannLevel
          cout << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ";
          cin >> config::RiemannLevel;
          break;
        case 'p': // Lev0Factor
          cout << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ";
          cin >> config::Lev0Factor;
          break;
        case 'r': // RTActualUse
          cout << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ";
          cin >> config::RTActualUse;
          break;
        case 's': // PsiType
          cout << Verbose(0) << "Old: " << config::PsiType << "\t new: ";
          cin >> config::PsiType;
          break;
        case 't': // MaxPsiDouble
          cout << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ";
          cin >> config::MaxPsiDouble;
          break;
        case 'u': // PsiMaxNoUp
          cout << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ";
          cin >> config::PsiMaxNoUp;
          break;
        case 'v': // PsiMaxNoDown
          cout << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ";
          cin >> config::PsiMaxNoDown;
          break;
        case 'w': // AddPsis
          cout << Verbose(0) << "Old: " << config::AddPsis << "\t new: ";
          cin >> config::AddPsis;
          break;

        case 'x': // RCut
          cout << Verbose(0) << "Old: " << config::RCut << "\t new: ";
          cin >> config::RCut;
          break;
        case 'y': // StructOpt
          cout << Verbose(0) << "Old: " << config::StructOpt << "\t new: ";
          cin >> config::StructOpt;
          break;
        case 'z': // IsAngstroem
          cout << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ";
          cin >> config::IsAngstroem;
          break;
        case 'i': // RelativeCoord
          cout << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ";
          cin >> config::RelativeCoord;
          break;
    };
  } while (choice != 'q');
};

/** Tests whether a given configuration file adhears to old or new syntax.
 * \param *filename filename of config file to be tested
 * \param *periode pointer to a periodentafel class with all elements
 * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax
 */
int config::TestSyntax(const char * const filename, const periodentafel * const periode) const
{
  int test;
  ifstream file(filename);

  // search file for keyword: ProcPEGamma (new syntax)
  if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
    file.close();
    return 1;
  }
  // search file for keyword: ProcsGammaPsi (old syntax)
  if (ParseForParameter(1,&file,"ProcsGammaPsi", 0, 1, 1, int_type, &test, 1, optional)) {
    file.close();
    return 0;
  }
  file.close();
  return -1;
}

/** Returns private config::IsAngstroem.
 * \return IsAngstroem
 */
bool config::GetIsAngstroem() const
{
  return (IsAngstroem == 1);
};

/** Returns private config::*defaultpath.
 * \return *defaultpath
 */
char * config::GetDefaultPath() const
{
  return defaultpath;
};


/** Returns private config::*defaultpath.
 * \return *defaultpath
 */
void config::SetDefaultPath(const char * const path)
{
  strcpy(defaultpath, path);
};

/** Retrieves the path in the given config file name.
 * \param filename config file string
 */
void config::RetrieveConfigPathAndName(const string filename)
{
  char *ptr = NULL;
  char *buffer = new char[MAXSTRINGSIZE];
  strncpy(buffer, filename.c_str(), MAXSTRINGSIZE);
  int last = -1;
  for(last=MAXSTRINGSIZE;last--;) {
    if (buffer[last] == '/')
      break;
  }
  if (last == -1) { // no path in front, set to local directory.
    strcpy(configpath, "./");
    ptr = buffer;
  } else {
    strncpy(configpath, buffer, last+1);
    ptr = &buffer[last+1];
    if (last < 254)
      configpath[last+1]='\0';
  }
  strcpy(configname, ptr);
  cout << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl;
  delete[](buffer);
};

/** Initializes ConfigFileBuffer from a file.
 * \param *file input file stream being the opened config file
 * \param *FileBuffer pointer to FileBuffer on return, should point to NULL
 */
void PrepareFileBuffer(const char * const filename, struct ConfigFileBuffer *&FileBuffer)
{
  if (FileBuffer != NULL) {
    cerr << Verbose(1) << "WARNING: deleting present FileBuffer in PrepareFileBuffer()." << endl;
    delete(FileBuffer);
  }
  FileBuffer = new ConfigFileBuffer(filename);

  FileBuffer->InitMapping();
};

/** Loads a molecule from a ConfigFileBuffer.
 * \param *mol molecule to load
 * \param *FileBuffer ConfigFileBuffer to use
 * \param *periode periodentafel for finding elements
 * \param FastParsing whether to parse trajectories or not
 */
void LoadMolecule(molecule * const &mol, struct ConfigFileBuffer * const &FileBuffer, const periodentafel * const periode, const bool FastParsing)
{
  int MaxTypes = 0;
  element *elementhash[MAX_ELEMENTS];
  char name[MAX_ELEMENTS];
  char keyword[MAX_ELEMENTS];
  int Z = -1;
  int No[MAX_ELEMENTS];
  int verbose = 0;
  double value[3];

  if (mol == NULL) {
    cerr << "Molecule is not allocated in LoadMolecule(), exit.";
    performCriticalExit();
  }

  ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
  if (MaxTypes == 0) {
    cerr << "There are no atoms according to MaxTypes in this config file." << endl;
  } else {
    // prescan number of ions per type
    cout << Verbose(0) << "Prescanning ions per type: " << endl;
    int NoAtoms = 0;
    for (int i=0; i < MaxTypes; i++) {
      sprintf(name,"Ion_Type%i",i+1);
      ParseForParameter(verbose,FileBuffer, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
      ParseForParameter(verbose,FileBuffer, name, 0, 2, 1, int_type, &Z, 1, critical);
      elementhash[i] = periode->FindElement(Z);
      cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
      NoAtoms += No[i];
    }
    int repetition = 0; // which repeated keyword shall be read

    // sort the lines via the LineMapping
    sprintf(name,"Ion_Type%i",MaxTypes);
    if (!ParseForParameter(verbose,FileBuffer, (const char*)name, 1, 1, 1, int_type, &value[0], 1, critical)) {
      cerr << "There are no atoms in the config file!" << endl;
      return;
    }
    FileBuffer->CurrentLine++;
    //cout << FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine]];
    FileBuffer->MapIonTypesInBuffer(NoAtoms);
    //for (int i=0; i<(NoAtoms < 100 ? NoAtoms : 100 < 100 ? NoAtoms : 100);++i) {
    //  cout << FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine+i]];
    //}

    map<int, atom *> AtomList[MaxTypes];
    map<int, atom *> LinearList;
    atom *neues = NULL;
    if (!FastParsing) {
      // parse in trajectories
      bool status = true;
      while (status) {
        cout << "Currently parsing MD step " << repetition << "." << endl;
        for (int i=0; i < MaxTypes; i++) {
          sprintf(name,"Ion_Type%i",i+1);
          for(int j=0;j<No[i];j++) {
            sprintf(keyword,"%s_%i",name, j+1);
            if (repetition == 0) {
              neues = new atom();
              AtomList[i][j] = neues;
              LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
              neues->type = elementhash[i]; // find element type
            } else
              neues = AtomList[i][j];
            status = (status &&
                    ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
                    ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
                    ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
                    ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
            if (!status) break;

            // check size of vectors
            if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
              //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
              neues->Trajectory.R.resize(repetition+10);
              neues->Trajectory.U.resize(repetition+10);
              neues->Trajectory.F.resize(repetition+10);
            }

            // put into trajectories list
            for (int d=0;d<NDIM;d++)
              neues->Trajectory.R.at(repetition).x[d] = neues->x.x[d];

            // parse velocities if present
            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
              neues->v.x[0] = 0.;
            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
              neues->v.x[1] = 0.;
            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
              neues->v.x[2] = 0.;
            for (int d=0;d<NDIM;d++)
              neues->Trajectory.U.at(repetition).x[d] = neues->v.x[d];

            // parse forces if present
            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
              value[0] = 0.;
            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
              value[1] = 0.;
            if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
              value[2] = 0.;
            for (int d=0;d<NDIM;d++)
              neues->Trajectory.F.at(repetition).x[d] = value[d];

  //            cout << "Parsed position of step " << (repetition) << ": (";
  //            for (int d=0;d<NDIM;d++)
  //              cout << neues->Trajectory.R.at(repetition).x[d] << " ";          // next step
  //            cout << ")\t(";
  //            for (int d=0;d<NDIM;d++)
  //              cout << neues->Trajectory.U.at(repetition).x[d] << " ";          // next step
  //            cout << ")\t(";
  //            for (int d=0;d<NDIM;d++)
  //              cout << neues->Trajectory.F.at(repetition).x[d] << " ";          // next step
  //            cout << ")" << endl;
          }
        }
        repetition++;
      }
      repetition--;
      cout << "Found " << repetition << " trajectory steps." << endl;
      if (repetition <= 1)  // if onyl one step, desactivate use of trajectories
        mol->MDSteps = 0;
      else
        mol->MDSteps = repetition;
    } else {
      // 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)
      repetition = 0;
      while ( ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
              ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
              ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
        repetition++;
      cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
      // parse in molecule coordinates
      for (int i=0; i < MaxTypes; i++) {
        sprintf(name,"Ion_Type%i",i+1);
        for(int j=0;j<No[i];j++) {
          sprintf(keyword,"%s_%i",name, j+1);
          if (repetition == 0) {
            neues = new atom();
            AtomList[i][j] = neues;
            LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
            neues->type = elementhash[i]; // find element type
          } else
            neues = AtomList[i][j];
          // then parse for each atom the coordinates as often as present
          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
          ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
            neues->v.x[0] = 0.;
          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
            neues->v.x[1] = 0.;
          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
            neues->v.x[2] = 0.;
          // here we don't care if forces are present (last in trajectories is always equal to current position)
          neues->type = elementhash[i]; // find element type
          mol->AddAtom(neues);
        }
      }
    }
    // put atoms into the molecule in their original order
    for(map<int, atom*>::iterator runner = LinearList.begin(); runner != LinearList.end(); ++runner) {
      mol->AddAtom(runner->second);
    }
  }
};


/** Initializes config file structure by loading elements from a give file.
 * \param *file input file stream being the opened config file
 * \param BondGraphFileName file name of the bond length table file, if string is left blank, no table is parsed.
 * \param *periode pointer to a periodentafel class with all elements
 * \param *&MolList pointer to MoleculeListClass, on return containing all parsed molecules in system
 */
void config::Load(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList)
{
  molecule *mol = new molecule(periode);
  ifstream *file = new ifstream(filename);
  if (file == NULL) {
    cerr << "ERROR: config file " << filename << " missing!" << endl;
    return;
  }
  file->close();
  delete(file);
  RetrieveConfigPathAndName(filename);

  // ParseParameterFile
  struct ConfigFileBuffer *FileBuffer = NULL;
  PrepareFileBuffer(filename,FileBuffer);

  /* Oeffne Hauptparameterdatei */
  int di = 0;
  double BoxLength[9];
  string zeile;
  string dummy;
  int verbose = 0;
  
  ParseThermostats(FileBuffer);
  
  /* Namen einlesen */

  // 1. parse in options
  ParseForParameter(verbose,FileBuffer, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
  ParseForParameter(verbose,FileBuffer, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
  ParseForParameter(verbose,FileBuffer, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
  ParseForParameter(verbose,FileBuffer,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
  ParseForParameter(verbose,FileBuffer,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical);

  if (!ParseForParameter(verbose,FileBuffer,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional))
    config::Seed = 1;

  if(!ParseForParameter(verbose,FileBuffer,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) {
    config::DoOutOrbitals = 0;
  } else {
    if (config::DoOutOrbitals < 0) config::DoOutOrbitals = 0;
    if (config::DoOutOrbitals > 1) config::DoOutOrbitals = 1;
  }
  ParseForParameter(verbose,FileBuffer,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
  if (config::DoOutVis < 0) config::DoOutVis = 0;
  if (config::DoOutVis > 1) config::DoOutVis = 1;
  if (!ParseForParameter(verbose,FileBuffer,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional))
    config::VectorPlane = -1;
  if (!ParseForParameter(verbose,FileBuffer,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional))
    config::VectorCut = 0.;
  ParseForParameter(verbose,FileBuffer,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
  if (config::DoOutMes < 0) config::DoOutMes = 0;
  if (config::DoOutMes > 1) config::DoOutMes = 1;
  if (!ParseForParameter(verbose,FileBuffer,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional))
    config::DoOutCurrent = 0;
  if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
  if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
  ParseForParameter(verbose,FileBuffer,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
  if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
  if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
  if(!ParseForParameter(verbose,FileBuffer,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) {
    config::DoWannier = 0;
  } else {
    if (config::DoWannier < 0) config::DoWannier = 0;
    if (config::DoWannier > 1) config::DoWannier = 1;
  }
  if(!ParseForParameter(verbose,FileBuffer,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) {
    config::CommonWannier = 0;
  } else {
    if (config::CommonWannier < 0) config::CommonWannier = 0;
    if (config::CommonWannier > 4) config::CommonWannier = 4;
  }
  if(!ParseForParameter(verbose,FileBuffer,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) {
    config::SawtoothStart = 0.01;
  } else {
    if (config::SawtoothStart < 0.) config::SawtoothStart = 0.;
    if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
  }
  
  if (ParseForParameter(verbose,FileBuffer,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional))
    if (config::DoConstrainedMD < 0) 
      config::DoConstrainedMD = 0;
  ParseForParameter(verbose,FileBuffer,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
  if (!ParseForParameter(verbose,FileBuffer,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
    config::Deltat = 1;
  ParseForParameter(verbose,FileBuffer,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
  ParseForParameter(verbose,FileBuffer,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
  ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
  //ParseForParameter(verbose,FileBuffer,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
  if (!ParseForParameter(verbose,FileBuffer,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
    config::EpsWannier = 1e-8;

  // stop conditions
  //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
  ParseForParameter(verbose,FileBuffer,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
  if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;

  ParseForParameter(verbose,FileBuffer,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
  ParseForParameter(verbose,FileBuffer,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
  ParseForParameter(verbose,FileBuffer,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
  ParseForParameter(verbose,FileBuffer,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical);
  ParseForParameter(verbose,FileBuffer,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical);
  if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
  if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
  if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;

  ParseForParameter(verbose,FileBuffer,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
  ParseForParameter(verbose,FileBuffer,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
  ParseForParameter(verbose,FileBuffer,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
  ParseForParameter(verbose,FileBuffer,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
  ParseForParameter(verbose,FileBuffer,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical);
  if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
  if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
  if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;

  // Unit cell and magnetic field
  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
  mol->cell_size[0] = BoxLength[0];
  mol->cell_size[1] = BoxLength[3];
  mol->cell_size[2] = BoxLength[4];
  mol->cell_size[3] = BoxLength[6];
  mol->cell_size[4] = BoxLength[7];
  mol->cell_size[5] = BoxLength[8];
  //if (1) fprintf(stderr,"\n");

  ParseForParameter(verbose,FileBuffer,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional);
  ParseForParameter(verbose,FileBuffer,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional);
  if (!ParseForParameter(verbose,FileBuffer,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional))
    config::DoFullCurrent = 0;
  if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
  if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
  if (config::DoOutNICS < 0) config::DoOutNICS = 0;
  if (config::DoOutNICS > 2) config::DoOutNICS = 2;
  if (config::DoPerturbation == 0) {
    config::DoFullCurrent = 0;
    config::DoOutNICS = 0;
  }

  ParseForParameter(verbose,FileBuffer,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
  ParseForParameter(verbose,FileBuffer,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
  ParseForParameter(verbose,FileBuffer,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
  if (config::Lev0Factor < 2) {
    config::Lev0Factor = 2;
  }
  ParseForParameter(verbose,FileBuffer,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
  if (di >= 0 && di < 2) {
    config::RiemannTensor = di;
  } else {
    fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    exit(1);
  }
  switch (config::RiemannTensor) {
    case 0: //UseNoRT
      if (config::MaxLevel < 2) {
        config::MaxLevel = 2;
      }
      config::LevRFactor = 2;
      config::RTActualUse = 0;
      break;
    case 1: // UseRT
      if (config::MaxLevel < 3) {
        config::MaxLevel = 3;
      }
      ParseForParameter(verbose,FileBuffer,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
      if (config::RiemannLevel < 2) {
        config::RiemannLevel = 2;
      }
      if (config::RiemannLevel > config::MaxLevel-1) {
        config::RiemannLevel = config::MaxLevel-1;
      }
      ParseForParameter(verbose,FileBuffer,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
      if (config::LevRFactor < 2) {
        config::LevRFactor = 2;
      }
      config::Lev0Factor = 2;
      config::RTActualUse = 2;
      break;
  }
  ParseForParameter(verbose,FileBuffer,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
  if (di >= 0 && di < 2) {
    config::PsiType = di;
  } else {
    fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    exit(1);
  }
  switch (config::PsiType) {
  case 0: // SpinDouble
    ParseForParameter(verbose,FileBuffer,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
    ParseForParameter(verbose,FileBuffer,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
    break;
  case 1: // SpinUpDown
    if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
    ParseForParameter(verbose,FileBuffer,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
    ParseForParameter(verbose,FileBuffer,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
    ParseForParameter(verbose,FileBuffer,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
    break;
  }

  // IonsInitRead

  ParseForParameter(verbose,FileBuffer,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
  ParseForParameter(verbose,FileBuffer,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
  ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
  if (!ParseForParameter(verbose,FileBuffer,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional))
    config::RelativeCoord = 0;
  if (!ParseForParameter(verbose,FileBuffer,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
    config::StructOpt = 0;

  // 2. parse the bond graph file if given
  BG = new BondGraph(IsAngstroem);
  if (BG->LoadBondLengthTable((ofstream *)&cout, BondGraphFileName)) {
    cout << Verbose(0) << "Bond length table loaded successfully." << endl;
  } else {
    cout << Verbose(0) << "Bond length table loading failed." << endl;
  }

  // 3. parse the molecule in
  LoadMolecule(mol, FileBuffer, periode, FastParsing);
  mol->ActiveFlag = true;
  MolList->insert(mol);

  // 4. dissect the molecule into connected subgraphs
  BG->ConstructBondGraph((ofstream *)&cout, mol);

  delete(FileBuffer);
};

/** Initializes config file structure by loading elements from a give file with old pcp syntax.
 * \param *file input file stream being the opened config file with old pcp syntax
 * \param BondGraphFileName file name of the bond length table file, if string is left blank, no table is parsed.
 * \param *periode pointer to a periodentafel class with all elements
 * \param *&MolList pointer to MoleculeListClass, on return containing all parsed molecules in system
 */
void config::LoadOld(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList)
{
  molecule *mol = new molecule(periode);
  ifstream *file = new ifstream(filename);
  if (file == NULL) {
    cerr << "ERROR: config file " << filename << " missing!" << endl;
    return;
  }
  RetrieveConfigPathAndName(filename);
  // ParseParameters

  /* Oeffne Hauptparameterdatei */
  int l = 0;
  int i = 0;
  int di = 0;
  double a = 0.;
  double b = 0.;
  double BoxLength[9];
  string zeile;
  string dummy;
  element *elementhash[128];
  int Z = -1;
  int No = -1;
  int AtomNo = -1;
  int found = 0;
  int verbose = 0;

  mol->ActiveFlag = true;
  MolList->insert(mol);
  /* Namen einlesen */

  ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
  ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
  ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
  ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
  ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 2, 1, int_type, &(config::ProcPEPsi), 1, critical);
  config::Seed = 1;
  config::DoOutOrbitals = 0;
  ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
  if (config::DoOutVis < 0) config::DoOutVis = 0;
  if (config::DoOutVis > 1) config::DoOutVis = 1;
  config::VectorPlane = -1;
  config::VectorCut = 0.;
  ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
  if (config::DoOutMes < 0) config::DoOutMes = 0;
  if (config::DoOutMes > 1) config::DoOutMes = 1;
  config::DoOutCurrent = 0;
  ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
  if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
  if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
  config::CommonWannier = 0;
  config::SawtoothStart = 0.01;

  ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, double_type, &(config::MaxOuterStep), 1, critical);
  ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
  ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
  ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
  ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
  ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
  config::EpsWannier = 1e-8;

  // stop conditions
  //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
  ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
  if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;

  ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
  ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
  ParseForParameter(verbose,file,"MaxMinStep", 0, 3, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
  ParseForParameter(verbose,file,"MaxMinStep", 0, 4, 1, int_type, &(config::MaxMinStopStep), 1, critical);
  if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
  if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
  config::MaxMinGapStopStep = 1;

  ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
  ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
  ParseForParameter(verbose,file,"MaxInitMinStep", 0, 3, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
  ParseForParameter(verbose,file,"MaxInitMinStep", 0, 4, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
  if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
  if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
  config::InitMaxMinGapStopStep = 1;

  ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
  mol->cell_size[0] = BoxLength[0];
  mol->cell_size[1] = BoxLength[3];
  mol->cell_size[2] = BoxLength[4];
  mol->cell_size[3] = BoxLength[6];
  mol->cell_size[4] = BoxLength[7];
  mol->cell_size[5] = BoxLength[8];
  if (1) fprintf(stderr,"\n");
  config::DoPerturbation = 0;
  config::DoFullCurrent = 0;

  ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
  ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
  ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
  if (config::Lev0Factor < 2) {
    config::Lev0Factor = 2;
  }
  ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
  if (di >= 0 && di < 2) {
    config::RiemannTensor = di;
  } else {
    fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    exit(1);
  }
  switch (config::RiemannTensor) {
    case 0: //UseNoRT
      if (config::MaxLevel < 2) {
        config::MaxLevel = 2;
      }
      config::LevRFactor = 2;
      config::RTActualUse = 0;
      break;
    case 1: // UseRT
      if (config::MaxLevel < 3) {
        config::MaxLevel = 3;
      }
      ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
      if (config::RiemannLevel < 2) {
        config::RiemannLevel = 2;
      }
      if (config::RiemannLevel > config::MaxLevel-1) {
        config::RiemannLevel = config::MaxLevel-1;
      }
      ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
      if (config::LevRFactor < 2) {
        config::LevRFactor = 2;
      }
      config::Lev0Factor = 2;
      config::RTActualUse = 2;
      break;
  }
  ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
  if (di >= 0 && di < 2) {
    config::PsiType = di;
  } else {
    fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    exit(1);
  }
  switch (config::PsiType) {
  case 0: // SpinDouble
    ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
    config::AddPsis = 0;
    break;
  case 1: // SpinUpDown
    if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
    ParseForParameter(verbose,file,"MaxPsiUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
    ParseForParameter(verbose,file,"MaxPsiDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
    config::AddPsis = 0;
    break;
  }

  // IonsInitRead

  ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
  ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
  config::RelativeCoord = 0;
  config::StructOpt = 0;


  // 2. parse the bond graph file if given
  BG = new BondGraph(IsAngstroem);
  if (BG->LoadBondLengthTable((ofstream *)&cout, BondGraphFileName)) {
    cout << Verbose(0) << "Bond length table loaded successfully." << endl;
  } else {
    cout << Verbose(0) << "Bond length table loading failed." << endl;
  }

  // Routine from builder.cpp

  for (i=MAX_ELEMENTS;i--;)
    elementhash[i] = NULL;
  cout << Verbose(0) << "Parsing Ions ..." << endl;
  No=0;
  found = 0;
  while (getline(*file,zeile,'\n')) {
    if (zeile.find("Ions_Data") == 0) {
      cout << Verbose(1) << "found Ions_Data...begin parsing" << endl;
      found ++;
    }
    if (found > 0) {
      if (zeile.find("Ions_Data") == 0)
        getline(*file,zeile,'\n'); // read next line and parse this one
      istringstream input(zeile);
      input >> AtomNo;  // number of atoms
      input >> Z;       // atomic number
      input >> a;
      input >> l;
      input >> l;
      input >> b;     // element mass
      elementhash[No] = periode->FindElement(Z);
      cout << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:"  << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl;
      for(i=0;i<AtomNo;i++) {
        if (!getline(*file,zeile,'\n')) {// parse on and on
          cout << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl;
          // return 1;
        } else {
          //cout << Verbose(2) << "Reading line: " << zeile << endl;
        }
        istringstream input2(zeile);
        atom *neues = new atom();
        input2 >> neues->x.x[0]; // x
        input2 >> neues->x.x[1]; // y
        input2 >> neues->x.x[2]; // z
        input2 >> l;
        neues->type = elementhash[No]; // find element type
        mol->AddAtom(neues);
      }
      No++;
    }
  }
  file->close();
  delete(file);
};

/** Stores all elements of config structure from which they can be re-read.
 * \param *filename name of file
 * \param *periode pointer to a periodentafel class with all elements
 * \param *mol pointer to molecule containing all atoms of the molecule
 */
bool config::Save(const char * const filename, const periodentafel * const periode, molecule * const mol) const
{
  bool result = true;
  // bring MaxTypes up to date
  mol->CountElements();
  ofstream * const output = new ofstream(filename, ios::out);
  if (output != NULL) {
    *output << "# ParallelCarParinello - main configuration file - created with molecuilder" << endl;
    *output << endl;
    *output << "mainname\t" << config::mainname << "\t# programm name (for runtime files)" << endl;
    *output << "defaultpath\t" << config::defaultpath << "\t# where to put files during runtime" << endl;
    *output << "pseudopotpath\t" << config::pseudopotpath << "\t# where to find pseudopotentials" << endl;
    *output << endl;
    *output << "ProcPEGamma\t" << config::ProcPEGamma << "\t# for parallel computing: share constants" << endl;
    *output << "ProcPEPsi\t" << config::ProcPEPsi << "\t# for parallel computing: share wave functions" << endl;
    *output << "DoOutVis\t" << config::DoOutVis << "\t# Output data for OpenDX" << endl;
    *output << "DoOutMes\t" << config::DoOutMes << "\t# Output data for measurements" << endl;
    *output << "DoOutOrbitals\t" << config::DoOutOrbitals << "\t# Output all Orbitals" << endl;
    *output << "DoOutCurr\t" << config::DoOutCurrent << "\t# Ouput current density for OpenDx" << endl;
    *output << "DoOutNICS\t" << config::DoOutNICS << "\t# Output Nucleus independent current shieldings" << endl;
    *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
    *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
    *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
    switch(Thermostat) {
      default:
      case None:
        break;
      case Woodcock:
        *output << ScaleTempStep;
        break;
      case Gaussian:
        *output << ScaleTempStep;
        break;
      case Langevin:
        *output << TempFrequency << "\t" << alpha;
        break;
      case Berendsen:
        *output << TempFrequency;
        break;
      case NoseHoover:
        *output << HooverMass;
        break;
    };
    *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
    *output << "VectorPlane\t" << config::VectorPlane << "\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot" << endl;
    *output << "VectorCut\t" << config::VectorCut << "\t# Cut plane axis value" << endl;
    *output << "AddGramSch\t" << config::UseAddGramSch << "\t# Additional GramSchmidtOrtogonalization to be safe" << endl;
    *output << "Seed\t\t" << config::Seed << "\t# initial value for random seed for Psi coefficients" << endl;
    *output << endl;
    *output << "MaxOuterStep\t" << config::MaxOuterStep << "\t# number of MolecularDynamics/Structure optimization steps" << endl;
    *output << "Deltat\t" << config::Deltat << "\t# time per MD step" << endl;
    *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
    *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
    *output << "TargetTemp\t" << config::TargetTemp << "\t# Target temperature" << endl;
    *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
    *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
    *output << endl;
    *output << "# Values specifying when to stop" << endl;
    *output << "MaxMinStep\t" << config::MaxMinStep << "\t# Maximum number of steps" << endl;
    *output << "RelEpsTotalE\t" << config::RelEpsTotalEnergy << "\t# relative change in total energy" << endl;
    *output << "RelEpsKineticE\t" << config::RelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
    *output << "MaxMinStopStep\t" << config::MaxMinStopStep << "\t# check every ..th steps" << endl;
    *output << "MaxMinGapStopStep\t" << config::MaxMinGapStopStep << "\t# check every ..th steps" << endl;
    *output << endl;
    *output << "# Values specifying when to stop for INIT, otherwise same as above" << endl;
    *output << "MaxInitMinStep\t" << config::MaxInitMinStep << "\t# Maximum number of steps" << endl;
    *output << "InitRelEpsTotalE\t" << config::InitRelEpsTotalEnergy << "\t# relative change in total energy" << endl;
    *output << "InitRelEpsKineticE\t" << config::InitRelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
    *output << "InitMaxMinStopStep\t" << config::InitMaxMinStopStep << "\t# check every ..th steps" << endl;
    *output << "InitMaxMinGapStopStep\t" << config::InitMaxMinGapStopStep << "\t# check every ..th steps" << endl;
    *output << endl;
    *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    *output << mol->cell_size[0] << "\t" << endl;
    *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
    *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
    // FIXME
    *output << endl;
    *output << "ECut\t\t" << config::ECut << "\t# energy cutoff for discretization in Hartrees" << endl;
    *output << "MaxLevel\t" << config::MaxLevel << "\t# number of different levels in the code, >=2" << endl;
    *output << "Level0Factor\t" << config::Lev0Factor << "\t# factor by which node number increases from S to 0 level" << endl;
    *output << "RiemannTensor\t" << config::RiemannTensor << "\t# (Use metric)" << endl;
    switch (config::RiemannTensor) {
      case 0: //UseNoRT
        break;
      case 1: // UseRT
        *output << "RiemannLevel\t" << config::RiemannLevel << "\t# Number of Riemann Levels" << endl;
        *output << "LevRFactor\t" << config::LevRFactor << "\t# factor by which node number increases from 0 to R level from" << endl;
        break;
    }
    *output << "PsiType\t\t" << config::PsiType << "\t# 0 - doubly occupied, 1 - SpinUp,SpinDown" << endl;
    // write out both types for easier changing afterwards
  //  switch (PsiType) {
  //    case 0:
        *output << "MaxPsiDouble\t" << config::MaxPsiDouble << "\t# here: specifying both maximum number of SpinUp- and -Down-states" << endl;
  //      break;
  //    case 1:
        *output << "PsiMaxNoUp\t" << config::PsiMaxNoUp << "\t# here: specifying maximum number of SpinUp-states" << endl;
        *output << "PsiMaxNoDown\t" << config::PsiMaxNoDown << "\t# here: specifying maximum number of SpinDown-states" << endl;
  //      break;
  //  }
    *output << "AddPsis\t\t" << config::AddPsis << "\t# Additional unoccupied Psis for bandgap determination" << endl;
    *output << endl;
    *output << "RCut\t\t" << config::RCut << "\t# R-cut for the ewald summation" << endl;
    *output << "StructOpt\t" << config::StructOpt << "\t# Do structure optimization beforehand" << endl;
    *output << "IsAngstroem\t" << config::IsAngstroem << "\t# 0 - Bohr, 1 - Angstroem" << endl;
    *output << "RelativeCoord\t" << config::RelativeCoord << "\t# whether ion coordinates are relative (1) or absolute (0)" << endl;
    *output << "MaxTypes\t" << mol->ElementCount <<  "\t# maximum number of different ion types" << endl;
    *output << endl;
    result = result && mol->Checkout(output);
    if (mol->MDSteps <=1 )
      result = result && mol->Output(output);
    else
      result = result && mol->OutputTrajectories(output);
    output->close();
    output->clear();
    delete(output);
    return result;
  } else
    return false;
};

/** Stores all elements in a MPQC input file.
 * Note that this format cannot be parsed again.
 * \param *filename name of file (without ".in" suffix!)
 * \param *mol pointer to molecule containing all atoms of the molecule
 */
bool config::SaveMPQC(const char * const filename, const molecule * const mol) const
{
  int AtomNo = -1;
  Vector *center = NULL;
  ofstream *output = NULL;

  // first without hessian
  {
    stringstream * const fname = new stringstream;;
    *fname << filename << ".in";
    output = new ofstream(fname->str().c_str(), ios::out);
    *output << "% Created by MoleCuilder" << endl;
    *output << "mpqc: (" << endl;
    *output << "\tsavestate = no" << endl;
    *output << "\tdo_gradient = yes" << endl;
    *output << "\tmole<MBPT2>: (" << endl;
    *output << "\t\tmaxiter = 200" << endl;
    *output << "\t\tbasis = $:basis" << endl;
    *output << "\t\tmolecule = $:molecule" << endl;
    *output << "\t\treference<CLHF>: (" << endl;
    *output << "\t\t\tbasis = $:basis" << endl;
    *output << "\t\t\tmolecule = $:molecule" << endl;
    *output << "\t\t)" << endl;
    *output << "\t)" << endl;
    *output << ")" << endl;
    *output << "molecule<Molecule>: (" << endl;
    *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
    *output << "\t{ atoms geometry } = {" << endl;
    center = mol->DetermineCenterOfAll(output);
    // output of atoms
    AtomNo = 0;
    mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
    delete(center);
    *output << "\t}" << endl;
    *output << ")" << endl;
    *output << "basis<GaussianBasisSet>: (" << endl;
    *output << "\tname = \"" << basis << "\"" << endl;
    *output << "\tmolecule = $:molecule" << endl;
    *output << ")" << endl;
    output->close();
    delete(output);
    delete(fname);
  }

  // second with hessian
  {
    stringstream * const fname = new stringstream;
    *fname << filename << ".hess.in";
    output = new ofstream(fname->str().c_str(), ios::out);
    *output << "% Created by MoleCuilder" << endl;
    *output << "mpqc: (" << endl;
    *output << "\tsavestate = no" << endl;
    *output << "\tdo_gradient = yes" << endl;
    *output << "\tmole<CLHF>: (" << endl;
    *output << "\t\tmaxiter = 200" << endl;
    *output << "\t\tbasis = $:basis" << endl;
    *output << "\t\tmolecule = $:molecule" << endl;
    *output << "\t)" << endl;
    *output << "\tfreq<MolecularFrequencies>: (" << endl;
    *output << "\t\tmolecule=$:molecule" << endl;
    *output << "\t)" << endl;
    *output << ")" << endl;
    *output << "molecule<Molecule>: (" << endl;
    *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
    *output << "\t{ atoms geometry } = {" << endl;
    center = mol->DetermineCenterOfAll(output);
    // output of atoms
    AtomNo = 0;
    mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
    delete(center);
    *output << "\t}" << endl;
    *output << ")" << endl;
    *output << "basis<GaussianBasisSet>: (" << endl;
    *output << "\tname = \"3-21G\"" << endl;
    *output << "\tmolecule = $:molecule" << endl;
    *output << ")" << endl;
    output->close();
    delete(output);
    delete(fname);
  }

  return true;
};

/** Reads parameter from a parsed file.
 * The file is either parsed for a certain keyword or if null is given for
 * the value in row yth and column xth. If the keyword was necessity#critical,
 * then an error is thrown and the programme aborted.
 * \warning value is modified (both in contents and position)!
 * \param verbose 1 - print found value to stderr, 0 - don't
 * \param *file file to be parsed
 * \param name Name of value in file (at least 3 chars!)
 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
 *        (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
 *         best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
 *         counted from this unresetted position!)
 * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
 * \param yth In grid case specifying column number, otherwise the yth \a name matching line
 * \param type Type of the Parameter to be read
 * \param value address of the value to be read (must have been allocated)
 * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
 * \param critical necessity of this keyword being specified (optional, critical)
 * \return 1 - found, 0 - not found
 * \note Routine is taken from the pcp project and hack-a-slack adapted to C++
 */
int ParseForParameter(const int verbose, ifstream * const file, const char * const name, const int sequential, const int xth, const int yth, const int type, void * value, const int repetition, const int critical) {
  int i = 0;
  int j = 0;  // loop variables
  int length = 0;
  int maxlength = -1;
  long file_position = file->tellg(); // mark current position
  char *dummy1 = NULL;
  char *dummy = NULL;
  char * const free_dummy = Malloc<char>(256, "config::ParseForParameter: *free_dummy");    // pointers in the line that is read in per step
  dummy1 = free_dummy;

  //fprintf(stderr,"Parsing for %s\n",name);
  if (repetition == 0)
    //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
    return 0;

  int line = 0; // marks line where parameter was found
  int found = (type >= grid) ? 0 : (-yth + 1);  // marks if yth parameter name was found
  while((found != repetition)) {
    dummy1 = dummy = free_dummy;
    do {
      file->getline(dummy1, 256); // Read the whole line
      if (file->eof()) {
        if ((critical) && (found == 0)) {
          Free(free_dummy);
          //Error(InitReading, name);
          fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
          exit(255);
        } else {
          //if (!sequential)
          file->clear();
          file->seekg(file_position, ios::beg);  // rewind to start position
          Free(free_dummy);
          return 0;
        }
      }
      line++;
    } while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines

    // C++ getline removes newline at end, thus re-add
    if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
      i = strlen(dummy1);
      dummy1[i] = '\n';
      dummy1[i+1] = '\0';
    }
    //fprintf(stderr,"line %i ends at %i, newline at %i\n", line, strlen(dummy1), strchr(dummy1,'\n')-free_dummy);

    if (dummy1 == NULL) {
      if (verbose) fprintf(stderr,"Error reading line %i\n",line);
    } else {
      //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
    }
    // Seek for possible end of keyword on line if given ...
    if (name != NULL) {
      dummy = strchr(dummy1,'\t');  // set dummy on first tab or space which ever's nearer
      if (dummy == NULL) {
        dummy = strchr(dummy1, ' ');  // if not found seek for space
        while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
          dummy++;
      }
      if (dummy == NULL) {
        dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
        //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
        //Free((void **)&free_dummy);
        //Error(FileOpenParams, NULL);
      } else {
        //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
      }
    } else dummy = dummy1;
    // ... and check if it is the keyword!
    //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
    if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
      found++; // found the parameter!
      //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);

      if (found == repetition) {
        for (i=0;i<xth;i++) { // i = rows
          if (type >= grid) {
            // grid structure means that grid starts on the next line, not right after keyword
            dummy1 = dummy = free_dummy;
            do {
              file->getline(dummy1, 256); // Read the whole line, skip commentary and empty ones
              if (file->eof()) {
                if ((critical) && (found == 0)) {
                  Free(free_dummy);
                  //Error(InitReading, name);
                  fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
                  exit(255);
                } else {
                  //if (!sequential)
                  file->clear();
                  file->seekg(file_position, ios::beg);  // rewind to start position
                  Free(free_dummy);
                  return 0;
                }
              }
              line++;
            } while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
            if (dummy1 == NULL){
              if (verbose) fprintf(stderr,"Error reading line %i\n", line);
            } else {
              //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
            }
          } else { // simple int, strings or doubles start in the same line
            while ((*dummy == '\t') || (*dummy == ' '))   // skip interjacent tabs and spaces
              dummy++;
          }
          // C++ getline removes newline at end, thus re-add
          if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
            j = strlen(dummy1);
            dummy1[j] = '\n';
            dummy1[j+1] = '\0';
          }

          int start = (type >= grid) ? 0 : yth-1 ;
          for (j=start;j<yth;j++) { // j = columns
            // check for lower triangular area and upper triangular area
            if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
              *((double *)value) = 0.0;
              fprintf(stderr,"%f\t",*((double *)value));
              value = (void *)((long)value + sizeof(double));
              //value += sizeof(double);
            } else {
              // otherwise we must skip all interjacent tabs and spaces and find next value
              dummy1 = dummy;
              dummy = strchr(dummy1, '\t'); // seek for tab or space
              if (dummy == NULL)
                dummy = strchr(dummy1, ' ');  // if not found seek for space
              if (dummy == NULL) { // if still zero returned ...
                dummy = strchr(dummy1, '\n');  // ... at line end then
                if ((j < yth-1) && (type < 4)) {  // check if xth value or not yet
                  if (critical) {
                    if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
                    Free(free_dummy);
                    //return 0;
                    exit(255);
                    //Error(FileOpenParams, NULL);
                  } else {
                    //if (!sequential)
                    file->clear();
                    file->seekg(file_position, ios::beg);  // rewind to start position
                    Free(free_dummy);
                    return 0;
                  }
                }
              } else {
                //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
              }
              if (*dummy1 == '#') {
                // found comment, skipping rest of line
                //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
                if (!sequential) { // here we need it!
                  file->seekg(file_position, ios::beg);  // rewind to start position
                }
                Free(free_dummy);
                return 0;
              }
              //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
              switch(type) {
                case (row_int):
                  *((int *)value) = atoi(dummy1);
                  if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
                  if (verbose) fprintf(stderr,"%i\t",*((int *)value));
                    value = (void *)((long)value + sizeof(int));
                    //value += sizeof(int);
                  break;
                case(row_double):
                case(grid):
                case(lower_trigrid):
                case(upper_trigrid):
                  *((double *)value) = atof(dummy1);
                  if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
                  if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
                  value = (void *)((long)value + sizeof(double));
                  //value += sizeof(double);
                  break;
                case(double_type):
                  *((double *)value) = atof(dummy1);
                  if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
                  //value += sizeof(double);
                  break;
                case(int_type):
                  *((int *)value) = atoi(dummy1);
                  if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
                  //value += sizeof(int);
                  break;
                default:
                case(string_type):
                  if (value != NULL) {
                    //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
                    maxlength = MAXSTRINGSIZE;
                    length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
                    strncpy((char *)value, dummy1, length);  // copy as much
                    ((char *)value)[length] = '\0';  // and set end marker
                    if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
                    //value += sizeof(char);
                  } else {
                  }
                break;
              }
            }
            while (*dummy == '\t')
              dummy++;
          }
        }
      }
    }
  }
  if ((type >= row_int) && (verbose))
    fprintf(stderr,"\n");
  Free(free_dummy);
  if (!sequential) {
    file->clear();
    file->seekg(file_position, ios::beg);  // rewind to start position
  }
  //fprintf(stderr, "End of Parsing\n\n");

  return (found); // true if found, false if not
}


/** Reads parameter from a parsed file.
 * The file is either parsed for a certain keyword or if null is given for
 * the value in row yth and column xth. If the keyword was necessity#critical,
 * then an error is thrown and the programme aborted.
 * \warning value is modified (both in contents and position)!
 * \param verbose 1 - print found value to stderr, 0 - don't
 * \param *FileBuffer pointer to buffer structure
 * \param name Name of value in file (at least 3 chars!)
 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
 *        (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
 *         best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
 *         counted from this unresetted position!)
 * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
 * \param yth In grid case specifying column number, otherwise the yth \a name matching line
 * \param type Type of the Parameter to be read
 * \param value address of the value to be read (must have been allocated)
 * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
 * \param critical necessity of this keyword being specified (optional, critical)
 * \return 1 - found, 0 - not found
 * \note Routine is taken from the pcp project and hack-a-slack adapted to C++
 */
int ParseForParameter(const int verbose, struct ConfigFileBuffer * const FileBuffer, const char * const name, const int sequential, const int xth, const int yth, const int type, void * value, const int repetition, const int critical) {
  int i = 0;
  int j = 0;  // loop variables
  int length = 0;
  int maxlength = -1;
  int OldCurrentLine = FileBuffer->CurrentLine;
  char *dummy1 = NULL;
  char *dummy = NULL;    // pointers in the line that is read in per step

  //fprintf(stderr,"Parsing for %s\n",name);
  if (repetition == 0)
    //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
    return 0;

  int line = 0; // marks line where parameter was found
  int found = (type >= grid) ? 0 : (-yth + 1);  // marks if yth parameter name was found
  while((found != repetition)) {
    dummy1 = dummy = NULL;
    do {
      dummy1 = FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine++] ];
      if (FileBuffer->CurrentLine >= FileBuffer->NoLines) {
        if ((critical) && (found == 0)) {
          //Error(InitReading, name);
          fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
          exit(255);
        } else {
          FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
          return 0;
        }
      }
      if (dummy1 == NULL) {
        if (verbose) fprintf(stderr,"Error reading line %i\n",line);
      } else {
        //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
      }
      line++;
    } while (dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines

    // Seek for possible end of keyword on line if given ...
    if (name != NULL) {
      dummy = strchr(dummy1,'\t');  // set dummy on first tab or space which ever's nearer
      if (dummy == NULL) {
        dummy = strchr(dummy1, ' ');  // if not found seek for space
        while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
          dummy++;
      }
      if (dummy == NULL) {
        dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
        //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
        //Free(&free_dummy);
        //Error(FileOpenParams, NULL);
      } else {
        //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
      }
    } else dummy = dummy1;
    // ... and check if it is the keyword!
    //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
    if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
      found++; // found the parameter!
      //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);

      if (found == repetition) {
        for (i=0;i<xth;i++) { // i = rows
          if (type >= grid) {
            // grid structure means that grid starts on the next line, not right after keyword
            dummy1 = dummy = NULL;
            do {
              dummy1 = FileBuffer->buffer[ FileBuffer->LineMapping[ FileBuffer->CurrentLine++] ];
              if (FileBuffer->CurrentLine >= FileBuffer->NoLines) {
                if ((critical) && (found == 0)) {
                  //Error(InitReading, name);
                  fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
                  exit(255);
                } else {
                  FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
                  return 0;
                }
              }
              if (dummy1 == NULL) {
                if (verbose) fprintf(stderr,"Error reading line %i\n", line);
              } else {
                //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
              }
              line++;
            } while (dummy1 != NULL && (dummy1[0] == '#') || (dummy1[0] == '\n'));
            dummy = dummy1;
          } else { // simple int, strings or doubles start in the same line
            while ((*dummy == '\t') || (*dummy == ' '))  // skip interjacent tabs and spaces
              dummy++;
          }

          for (j=((type >= grid) ? 0 : yth-1);j<yth;j++) { // j = columns
            // check for lower triangular area and upper triangular area
            if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
              *((double *)value) = 0.0;
              fprintf(stderr,"%f\t",*((double *)value));
              value = (void *)((long)value + sizeof(double));
              //value += sizeof(double);
            } else {
              // otherwise we must skip all interjacent tabs and spaces and find next value
              dummy1 = dummy;
              dummy = strchr(dummy1, '\t'); // seek for tab or space
              if (dummy == NULL)
                dummy = strchr(dummy1, ' ');  // if not found seek for space
              if (dummy == NULL) { // if still zero returned ...
                dummy = strchr(dummy1, '\n'); // ... at line end then
                if ((j < yth-1) && (type < 4)) {  // check if xth value or not yet
                  if (critical) {
                    if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
                    //return 0;
                    exit(255);
                    //Error(FileOpenParams, NULL);
                  } else {
                    if (!sequential) { // here we need it!
                      FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
                    }
                    return 0;
                  }
                }
              } else {
                //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
              }
              if (*dummy1 == '#') {
                // found comment, skipping rest of line
                //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
                if (!sequential) { // here we need it!
                  FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
                }
                return 0;
              }
              //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
              switch(type) {
                case (row_int):
                  *((int *)value) = atoi(dummy1);
                  if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
                  if (verbose) fprintf(stderr,"%i\t",*((int *)value));
                    value = (void *)((long)value + sizeof(int));
                    //value += sizeof(int);
                  break;
                case(row_double):
                case(grid):
                case(lower_trigrid):
                case(upper_trigrid):
                  *((double *)value) = atof(dummy1);
                  if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
                  if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
                  value = (void *)((long)value + sizeof(double));
                  //value += sizeof(double);
                  break;
                case(double_type):
                  *((double *)value) = atof(dummy1);
                  if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
                  //value += sizeof(double);
                  break;
                case(int_type):
                  *((int *)value) = atoi(dummy1);
                  if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
                  //value += sizeof(int);
                  break;
                default:
                case(string_type):
                  if (value != NULL) {
                    //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
                    maxlength = MAXSTRINGSIZE;
                    length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
                    strncpy((char *)value, dummy1, length); // copy as much
                    ((char *)value)[length] = '\0'; // and set end marker
                    if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
                    //value += sizeof(char);
                  } else {
                  }
                break;
              }
            }
            while (*dummy == '\t')
              dummy++;
          }
        }
      }
    }
  }
  if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
  if (!sequential) {
    FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
  }
  //fprintf(stderr, "End of Parsing\n\n");

  return (found); // true if found, false if not
}