source: src/config.cpp@ 47d041

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

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

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include <stdio.h>
#include <cstring>
#include <iostream>
#include <boost/tokenizer.hpp>
#include <string>

#include "atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/toString.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "ConfigFileBuffer.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "Graph/BondGraph.hpp"
#include "Helpers/helpers.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "ThermoStatContainer.hpp"
#include "World.hpp"

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

/** Constructor for config file class.
 */
config::config() :
  PsiType(0),
  MaxPsiDouble(0),
  PsiMaxNoUp(0),
  PsiMaxNoDown(0),
  MaxMinStopStep(1),
  InitMaxMinStopStep(1),
  ProcPEGamma(8),
  ProcPEPsi(1),
  configname(NULL),
  FastParsing(false),
  Deltat(0.01),
  databasepath(NULL),
  DoConstrainedMD(0),
  MaxOuterStep(0),
  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 = new char[MAXSTRINGSIZE];
  defaultpath = new char[MAXSTRINGSIZE];
  pseudopotpath = new char[MAXSTRINGSIZE];
  databasepath = new char[MAXSTRINGSIZE];
  configname = new char[MAXSTRINGSIZE];
  strcpy(mainname,"pcp");
  strcpy(defaultpath,"not specified");
  strcpy(pseudopotpath,"not specified");
  configname[0]='\0';
};

/** Destructor for config file class.
 */
config::~config()
{
  delete[](mainname);
  delete[](defaultpath);
  delete[](pseudopotpath);
  delete[](databasepath);
  delete[](configname);
};

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

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

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

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

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

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

//        case 'j': // BoxLength
//          std::cout << "enter lower triadiagonalo form of basis matrix" << std::endl << std::endl;
//          double * const cell_size = World::getInstance().getDomain();
//          for (int i=0;i<6;i++) {
//            std::cout << "Cell size" << i << ": ";
//            cin >> cell_size[i];
//          }
//          break;

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

        case 'x': // RCut
          std::cout << "Old: " << config::RCut << "\t new: " << std::endl;
          cin >> config::RCut;
          break;
        case 'y': // StructOpt
          std::cout << "Old: " << config::StructOpt << "\t new: " << std::endl;
          cin >> config::StructOpt;
          break;
        case 'z': // IsAngstroem
          std::cout << "Old: " << config::IsAngstroem << "\t new: " << std::endl;
          cin >> config::IsAngstroem;
          break;
        case 'i': // RelativeCoord
          std::cout << "Old: " << config::RelativeCoord << "\t new: " << std::endl;
          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);
};

/** 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;
  const element *elementhash[MAX_ELEMENTS];
  char name[MAXSTRINGSIZE];
  int Z = -1;
  int No[MAX_ELEMENTS];
  int verbose = DoLog(4);
  double value[3];

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

  ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
  if (MaxTypes == 0) {
    ELOG(1, "There are no atoms according to MaxTypes in this config file." << endl);
    //performCriticalExit();
  } else {
    // prescan number of ions per type
    LOG(0, "STATUS: 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);
      LOG(1, i << ". Z = " << elementhash[i]->getAtomicNumber() << " with " << No[i] << " ions.");
      NoAtoms += No[i];
    }
    int repetition = -1; // 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)) {
      ELOG(0, "There are no atoms in the config file!" << endl);
      performCriticalExit();
      return;
    }

    FileBuffer->CurrentLine++;  // skip to next line
    FileBuffer->MapIonTypesInBuffer(NoAtoms);
    for (int i=FileBuffer->CurrentLine; i<FileBuffer->NoLines;++i) {
      LOG(4, FileBuffer->buffer[ FileBuffer->LineMapping[i] ]);
    }

    map<int, atom *> AtomList[MaxTypes];
    map<int, atom *> LinearList;
    atom *neues = NULL;
    Vector tempVector;
    int _fixedion;

    typedef boost::tokenizer<boost::char_separator<char> >
        tokenizer;
    boost::char_separator<char> sep("\t ");
    ConvertTo<double> toDouble;
    ConvertTo<int> toInt;

    for (int i=0; i < MaxTypes; i++) {
      for(int j=0;j<No[i];j++) {
        int step = 0;
        std::stringstream keyword_stream;
        keyword_stream << "Ion_Type" << i+1 << "_" << j+1;
        const std::string keyword = keyword_stream.str();
        LOG(3, "INFO: Parsing for " << keyword << "." << std::endl);
        while (true) {
          const std::string line(FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ]);
          const std::string line_without_comment = line.substr(0,line.find("#"));
          tokenizer tokens(line_without_comment, sep);
          if (tokens.begin() != tokens.end()) {
            tokenizer::iterator tok_iter = tokens.begin();
            const std::string token = *tok_iter++;
            if (token == keyword) {
              LOG(3, "INFO: Found keyword " << keyword << " in line " << FileBuffer->CurrentLine << std::endl);
              if (step == 0) {
                neues = World::getInstance().createAtom();
                AtomList[i][j] = neues;
                LOG(4, "Filling LinearList [ (FileBuffer->LineMapping[" << FileBuffer->CurrentLine << "]) = " << FileBuffer->LineMapping[FileBuffer->CurrentLine] << " with " << neues << endl);
                LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
                neues->setType(elementhash[i]); // find element type
              } else
                neues = AtomList[i][j];

              // count tokens
              size_t tokens_size = 0;
              for (tokenizer::iterator tokiter = tokens.begin(); tokiter != tokens.end(); ++tokiter)
                ++tokens_size;
              LOG(3, "INFO: Line contains " << tokens_size << " tokens." << std::endl);
              // and parse
              tempVector.Zero();
              if (tokens_size >= 5) { // only AtomicPosition and FixedIon
                for (int i=0;i<NDIM;++i)
                  tempVector[i] = toDouble(*tok_iter++);
                neues->setPositionAtStep(step, tempVector);
                _fixedion = toInt(*tok_iter++);
                neues->setFixedIon(_fixedion == 1);
                LOG(3, "INFO: Parsing AtomicPosition " << tempVector << " and FixedIon " << _fixedion << "." << std::endl);
              }
              tempVector.Zero();
              if (tokens_size >= 8) { // AtomicVelocity
                for (int i=0;i<NDIM;++i)
                  tempVector[i] = toDouble(*tok_iter++);
                LOG(3, "INFO: Parsing AtomicVelocity " << tempVector << "." << std::endl);
              }
              neues->setAtomicVelocityAtStep(step, tempVector);
              tempVector.Zero();
              if (tokens_size >= 11) {  // AtomicForce
                LOG(3, "INFO: Parsing AtomicForce" << std::endl);
                for (int i=0;i<NDIM;++i)
                  tempVector[i] = toDouble(*tok_iter++);
              }
              neues->setAtomicForceAtStep(step, tempVector);
              std::stringstream output;
              output << "Parsed position of step " << (step+1) << ": ";
              output << neues->getPositionAtStep(step);          // next step
              output << "\t";
              output << (neues->getFixedIon() ? "true" : "false");
              output << "\t";
              output << neues->getAtomicVelocityAtStep(step);          // next step
              output << "\t";
              output << neues->getAtomicForceAtStep(step);          // next step
              LOG(2, output.str());

              step++;
            } else {
              if ((repetition > step) || (repetition == -1))
                repetition = step;
              break;
            }
          }
          FileBuffer->CurrentLine++;
        }
      }
    }

    if (repetition <= 1)  // if onyl one step, desactivate use of trajectories
      mol->MDSteps = 0;
    else {
      LOG(0, "Found " << repetition << " trajectory step(s).");
      mol->MDSteps = repetition;
    }

    // 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);
    }
  }
};

/** 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;
  const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
  ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
  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" << Thermostats->activeThermostat->name() << "\t";
    *output << Thermostats->activeThermostat->writeParams();
    *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" << Thermostats->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 << domain.at(0,0) << "\t" << endl;
    *output << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
    *output << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\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->getElementCount() <<  "\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 {
    ELOG(1, "Cannot open output file:" << filename);
    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
{
  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);
    if (output == NULL) {
      ELOG(1, "Cannot open mpqc output file:" << fname);
      delete(fname);
      return false;
    }
    *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 of atoms
    for(molecule::const_iterator iter = mol->begin(); iter!=mol->end();++iter){
      (*iter)->OutputMPQCLine(output,center);
    }
    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);
  }

  // second with hessian
  {
    stringstream * const fname = new stringstream;
    *fname << filename << ".hess.in";
    output = new ofstream(fname->str().c_str(), ios::out);
    if (output == NULL) {
      ELOG(1, "Cannot open mpqc hessian output file:" << fname);
      delete(fname);
      return false;
    }
    *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 of atoms
    for(molecule::const_iterator iter = mol->begin(); iter!=mol->end();++iter){
      (*iter)->OutputMPQCLine(output,center);
    }
    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;
};

/** Stores all atoms from all molecules in a PDB input file.
 * Note that this format cannot be parsed again.
 * \param *filename name of file (without ".in" suffix!)
 * \param *MolList pointer to MoleculeListClass containing all atoms
 */
bool config::SavePDB(const char * const filename, const MoleculeListClass * const MolList) const
{
  int AtomNo = -1;
  int MolNo = 0;
  FILE *f = NULL;

  char name[MAXSTRINGSIZE];
  strncpy(name, filename, MAXSTRINGSIZE-1);
  strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
  f = fopen(name, "w" );
  if (f == NULL) {
    ELOG(1, "Cannot open pdb output file:" << name);
    return false;
  }
  fprintf(f, "# Created by MoleCuilder\n");

  for (MoleculeList::const_iterator MolRunner = MolList->ListOfMolecules.begin(); MolRunner != MolList->ListOfMolecules.end(); MolRunner++) {
    int *elementNo = new int[MAX_ELEMENTS];
    for (int i=0;i<MAX_ELEMENTS;i++)
      elementNo[i] = 0;
    AtomNo = 0;
    for (molecule::const_iterator iter = (*MolRunner)->begin(); iter != (*MolRunner)->end(); ++iter) {
      sprintf(name, "%2s%2d",(*iter)->getType()->getSymbol().c_str(), elementNo[(*iter)->getType()->getAtomicNumber()]);
      elementNo[(*iter)->getType()->getAtomicNumber()] = (elementNo[(*iter)->getType()->getAtomicNumber()]+1) % 100;   // confine to two digits
      fprintf(f,
             "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
             (*iter)->getNr(),                /* atom serial number */
             name,         /* atom name */
             (*MolRunner)->name,      /* residue name */
             'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
             MolNo,         /* residue sequence number */
             (*iter)->at(0),                 /* position X in Angstroem */
             (*iter)->at(1),                 /* position Y in Angstroem */
             (*iter)->at(2),                 /* position Z in Angstroem */
             (double)(*iter)->getType()->getValence(),         /* occupancy */
             (double)(*iter)->getType()->getNoValenceOrbitals(),          /* temperature factor */
             "0",            /* segment identifier */
             (*iter)->getType()->getSymbol().c_str(),    /* element symbol */
             "0");           /* charge */
      AtomNo++;
    }
    delete[](elementNo);
    MolNo++;
  }
  fclose(f);

  return true;
};

/** Stores all atoms in a PDB input file.
 * Note that this format cannot be parsed again.
 * \param *filename name of file (without ".in" suffix!)
 * \param *mol pointer to molecule
 */
bool config::SavePDB(const char * const filename, const molecule * const mol) const
{
  int AtomNo = -1;
  FILE *f = NULL;

  int *elementNo = new int[MAX_ELEMENTS];
  for (int i=0;i<MAX_ELEMENTS;i++)
    elementNo[i] = 0;
  char name[MAXSTRINGSIZE];
  strncpy(name, filename, MAXSTRINGSIZE-1);
  strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
  f = fopen(name, "w" );
  if (f == NULL) {
    ELOG(1, "Cannot open pdb output file:" << name);
    delete[](elementNo);
    return false;
  }
  fprintf(f, "# Created by MoleCuilder\n");

  AtomNo = 0;
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    sprintf(name, "%2s%2d",(*iter)->getType()->getSymbol().c_str(), elementNo[(*iter)->getType()->getAtomicNumber()]);
    elementNo[(*iter)->getType()->getAtomicNumber()] = (elementNo[(*iter)->getType()->getAtomicNumber()]+1) % 100;   // confine to two digits
    fprintf(f,
           "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
           (*iter)->getNr(),                /* atom serial number */
           name,         /* atom name */
           mol->name,      /* residue name */
           'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
           0,         /* residue sequence number */
           (*iter)->at(0),                 /* position X in Angstroem */
           (*iter)->at(1),                 /* position Y in Angstroem */
           (*iter)->at(2),                 /* position Z in Angstroem */
           (double)(*iter)->getType()->getValence(),         /* occupancy */
           (double)(*iter)->getType()->getNoValenceOrbitals(),          /* temperature factor */
           "0",            /* segment identifier */
           (*iter)->getType()->getSymbol().c_str(),    /* element symbol */
           "0");           /* charge */
    AtomNo++;
  }
  fclose(f);
  delete[](elementNo);

  return true;
};

/** Stores all atoms in a TREMOLO data input file.
 * Note that this format cannot be parsed again.
 * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
 * \param *filename name of file (without ".in" suffix!)
 * \param *mol pointer to molecule
 */
bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
{
  ofstream *output = NULL;
  stringstream * const fname = new stringstream;

  *fname << filename << ".data";
  output = new ofstream(fname->str().c_str(), ios::out);
  if (output == NULL) {
    ELOG(1, "Cannot open tremolo output file:" << fname);
    delete(fname);
    return false;
  }

  // scan maximum number of neighbours
  int MaxNeighbours = 0;
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    const int count = (*iter)->getListOfBonds().size();
    if (MaxNeighbours < count)
      MaxNeighbours = count;
  }
  *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;

  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    *output << (*iter)->getNr() << "\t";
    *output << (*iter)->getName() << "\t";
    *output << mol->name << "\t";
    *output << 0 << "\t";
    *output << (*iter)->at(0) << "\t" << (*iter)->at(1) << "\t" << (*iter)->at(2) << "\t";
    *output << static_cast<double>((*iter)->getType()->getValence()) << "\t";
    *output << (*iter)->getType()->getSymbol() << "\t";
    const BondList& ListOfBonds = (*iter)->getListOfBonds();
    for (BondList::const_iterator runner = ListOfBonds.begin();
        runner != ListOfBonds.end();
        runner++) {
      *output << (*runner)->GetOtherAtom(*iter)->getNr() << "\t";
    }
    for(int i = ListOfBonds.size(); i < MaxNeighbours; i++)
      *output << "-\t";
    *output << endl;
  }
  output->flush();
  output->close();
  delete(output);
  delete(fname);

  return true;
};

/** Stores all atoms from all molecules in a TREMOLO data input file.
 * Note that this format cannot be parsed again.
 * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
 * \param *filename name of file (without ".in" suffix!)
 * \param *MolList pointer to MoleculeListClass containing all atoms
 */
bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
{
  Info FunctionInfo(__func__);
  ofstream *output = NULL;
  stringstream * const fname = new stringstream;

  *fname << filename << ".data";
  output = new ofstream(fname->str().c_str(), ios::out);
  if (output == NULL) {
    ELOG(1, "Cannot open tremolo output file:" << fname);
    delete(fname);
    return false;
  }

  // scan maximum number of neighbours
  int MaxNeighbours = 0;
  for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      const int count = (*iter)->getListOfBonds().size();
      if (MaxNeighbours < count)
        MaxNeighbours = count;
    }
  }
  *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;

  // create global to local id map
  map<int, int> LocalNotoGlobalNoMap;
  {
    unsigned int MolCounter = 0;
    int AtomNo = 1;
    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
      for(molecule::iterator AtomRunner = (*MolWalker)->begin(); AtomRunner != (*MolWalker)->end(); ++AtomRunner) {
        LocalNotoGlobalNoMap.insert( pair<int,int>((*AtomRunner)->getId(), AtomNo++) );
      }
      MolCounter++;
    }
    ASSERT(MolCounter == MolList->ListOfMolecules.size(), "SaveTREMOLO: LocalNotoGlobalNoMap[] has not been correctly initialized for each molecule");
  }

  // write the file
  {
    int MolCounter = 0;
    int AtomNo = 0;
    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        *output << LocalNotoGlobalNoMap[ (*iter)->getId() ] << "\t";
        *output << (*iter)->getName() << "\t";
        *output << (*MolWalker)->name << "\t";
        *output << MolCounter+1 << "\t";
        *output << (*iter)->at(0) << "\t" << (*iter)->at(1) << "\t" << (*iter)->at(2) << "\t";
        *output << (double)(*iter)->getType()->getValence() << "\t";
        *output << (*iter)->getType()->getSymbol() << "\t";
        const BondList& ListOfBonds = (*iter)->getListOfBonds();
        for (BondList::const_iterator runner = ListOfBonds.begin();
            runner != ListOfBonds.end();
            runner++) {
          *output << LocalNotoGlobalNoMap[ (*runner)->GetOtherAtom((*iter))->getId() ] << "\t";
        }
        for(int i = ListOfBonds.size(); i < MaxNeighbours; i++)
          *output << "-\t";
        *output << endl;
        AtomNo++;
      }
      MolCounter++;
    }
  }

  // store & free
  output->flush();
  output->close();
  delete(output);
  delete(fname);

  return true;
};


/** Tries given filename or standard on saving the config file.
 * \param *ConfigFileName name of file
 * \param *periode pointer to periodentafel structure with all the elements
 * \param *molecules list of molecules structure with all the atoms and coordinates
 */
void config::SaveAll(char *ConfigFileName, periodentafel *periode, MoleculeListClass *molecules)
{
  char filename[MAXSTRINGSIZE];
  ofstream output;
  molecule *mol = NULL;

  // first save as PDB data
  if (ConfigFileName != NULL)
    strcpy(filename, ConfigFileName);
  if (output == NULL)
    strcpy(filename,"main_pcp_linux");
  if (SavePDB(filename, molecules))
    LOG(0, "Saving as pdb input done.");
  else
    LOG(0, "Saving as pdb input failed.");

  // then save as tremolo data file
  if (ConfigFileName != NULL)
    strcpy(filename, ConfigFileName);
  if (output == NULL)
    strcpy(filename,"main_pcp_linux");
  if (SaveTREMOLO(filename, molecules))
    LOG(0, "Saving as tremolo data input done.");
  else
    LOG(0, "Saving as tremolo data input failed.");

  // translate each to its center and merge all molecules in MoleculeListClass into this molecule
  int N = molecules->ListOfMolecules.size();
  if (N != 1) { // don't do anything in case of only one molecule (shifts mol ids otherwise)
    int *src = new int[N];
    N=0;
    for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
      src[N++] = (*ListRunner)->IndexNr;
    }
    mol = World::getInstance().createMolecule();
    mol->SetNameFromFilename(ConfigFileName);
    //mol->CalculateOrbitals(*this);
    delete[](src);
  } else {
    if (!molecules->ListOfMolecules.empty()) {
      mol = *(molecules->ListOfMolecules.begin());
      mol->doCountAtoms();
      //mol->CalculateOrbitals(*this);
    } else {
      ELOG(1, "There are no molecules to save!");
    }
  }

  LOG(0, "Storing configuration ... ");
  // get correct valence orbitals
  if (ConfigFileName != NULL) { // test the file name
    strcpy(filename, ConfigFileName);
    output.open(filename, ios::trunc);
  } else if (strlen(configname) != 0) {
    strcpy(filename, configname);
    output.open(configname, ios::trunc);
    } else {
      strcpy(filename, DEFAULTCONFIG);
      output.open(DEFAULTCONFIG, ios::trunc);
    }
  output.close();
  output.clear();
  if (Save(filename, periode, mol))
    LOG(0, "Saving of config file done.");
  else
    LOG(0, "Saving of config file failed.");

  // and save to xyz file
  if (ConfigFileName != NULL) {
    strcpy(filename, ConfigFileName);
    strcat(filename, ".xyz");
    output.open(filename, ios::trunc);
  }
  if (output == NULL) {
    strcpy(filename,"main_pcp_linux");
    strcat(filename, ".xyz");
    output.open(filename, ios::trunc);
  }
  if (mol->MDSteps <= 1) {
    if (mol->OutputXYZ(&output))
      LOG(0, "Saving of XYZ file done.");
    else
      LOG(0, "Saving of XYZ file failed.");
  } else {
    if (mol->OutputTrajectoriesXYZ(&output))
      LOG(0, "Saving of XYZ file done.");
    else
      LOG(0, "Saving of XYZ file failed.");
  }
  output.close();
  output.clear();

  // and save as MPQC configuration
  if (ConfigFileName != NULL)
    strcpy(filename, ConfigFileName);
  if (output == NULL)
    strcpy(filename,"main_pcp_linux");
  if (SaveMPQC(filename, mol))
    LOG(0, "Saving as mpqc input done.");
  else
    LOG(0, "Saving as mpqc input failed.");

  // don't destroy molecule as it contains all our atoms
  //World::getInstance().destroyMolecule(mol);
};

/** 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 free_dummy[MAXSTRINGSIZE];    // 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)) {
          //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
          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);
        //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)) {
                  //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
                  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);
                    //return 0;
                    exit(255);
                    //Error(FileOpenParams, NULL);
                  } else {
                    //if (!sequential)
                    file->clear();
                    file->seekg(file_position, ios::beg);  // 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!
                  file->seekg(file_position, ios::beg);  // 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) {
    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
  char *free_dummy = NULL;

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

  int found = (type >= grid) ? 0 : (-yth + 1);  // marks if yth parameter name was found
  while((found != repetition)) {
    dummy1 = dummy = NULL;
    do {
      if (FileBuffer->CurrentLine < FileBuffer->NoLines)
        free_dummy = 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);
          return 0;
        } else {
          //fprintf(stdout,"Rewinding to OldCurrentLine due to search till end of file.\n");
          FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
          return 0;
        }
      }
      if (dummy1 == NULL) {
        if (verbose) fprintf(stderr,"Error reading line %i\n",FileBuffer->CurrentLine);
      } else {
        if (verbose) fprintf(stderr,"Now parsing the line %i: %s\n", FileBuffer->CurrentLine, dummy1);
      }
      //FileBuffer->CurrentLine++;
    } 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);
        //Error(FileOpenParams, NULL);
      } else {
        if (verbose) fprintf(stderr,"found tab at line %i at position %li\n",FileBuffer->CurrentLine, (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!
      if (verbose) fprintf(stderr,"found %s at line %i between %li and %li\n", name, FileBuffer->CurrentLine, (unsigned long)dummy1, (unsigned long)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 {
                  //fprintf(stdout,"Rewinding to OldCurrentLine due to search till end of line.\n");
                  FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
                  return 0;
                }
              }
              if (dummy1 == NULL) {
                if (verbose) fprintf(stderr,"Error reading line %i\n", FileBuffer->CurrentLine);
              } else {
                if (verbose) fprintf(stderr,"Reading next line %i: %s\n", FileBuffer->CurrentLine, dummy1);
              }
              //FileBuffer->CurrentLine++;
            } 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", FileBuffer->CurrentLine, yth-j, name);
                    //return 0;
                    exit(255);
                    //Error(FileOpenParams, NULL);
                  } else {
                    if (!sequential) { // here we need it!
                      //fprintf(stdout,"Rewinding to OldCurrentLine due to end of line and sequential %d.\n", sequential);
                      FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
                    }
                    return 0;
                  }
                }
              } else {
                if (verbose) fprintf(stderr,"found tab at line %i at position %li\n",FileBuffer->CurrentLine, (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!
                  //fprintf(stdout,"Rewinding to OldCurrentLine due to comment and sequential %d.\n", sequential);
                  FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
                }
                return 0;
              }
              if (verbose) fprintf(stderr,"value from %li to %li\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) {
    //fprintf(stdout,"Rewinding to OldCurrentLine due to sequential %d.\n", sequential);
    FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
  }
  if (verbose) fprintf(stderr, "End of Parsing for %s\n\n",name);

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

/** Reading 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 = new char[12];
  const int verbose = 0;
  ThermoStatContainer *Thermostats = World::getInstance().getThermostats();

  // read desired Thermostat from file along with needed additional parameters
  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    Thermostats->makeActive(thermo,fb);
  } else {
    if ((Thermostats->TargetTemp != 0))
      LOG(2,  "No thermostat chosen despite finite temperature MD, falling back to None.");
    Thermostats->chooseNone();
  }
  delete[](thermo);
};