/*
 * PcpParser.cpp
 *
 *  Created on: 12.06.2010
 *      Author: heber
 */

#include <iostream>

#include "atom.hpp"
#include "config.hpp"
#include "ConfigFileBuffer.hpp"
#include "element.hpp"
#include "log.hpp"
#include "molecule.hpp"
#include "PcpParser.hpp"
#include "periodentafel.hpp"
#include "ThermoStatContainer.hpp"
#include "verbose.hpp"
#include "World.hpp"

/** Constructor of PcpParser.
 *
 */
PcpParser::PcpParser()
{}

/** Destructor of PcpParser.
 *
 */
PcpParser::~PcpParser()
{}

void PcpParser::load(std::istream* file)
{
  if (file->fail()) {
    DoeLog(1) && (eLog()<< Verbose(1) << "could not access given file" << endl);
    return;
  }

  // ParseParameterFile
  class ConfigFileBuffer *FileBuffer = new ConfigFileBuffer();
  FileBuffer->InitFileBuffer(file);

  /* 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, (Paths.mainname), 1, critical);
  ParseForParameter(verbose,FileBuffer, "defaultpath", 0, 1, 1, string_type, (Paths.defaultpath), 1, critical);
  ParseForParameter(verbose,FileBuffer, "pseudopotpath", 0, 1, 1, string_type, (Paths.pseudopotpath), 1, critical);
  ParseForParameter(verbose,FileBuffer,"ProcPEGamma", 0, 1, 1, int_type, &(Parallelization.ProcPEGamma), 1, critical);
  ParseForParameter(verbose,FileBuffer,"ProcPEPsi", 0, 1, 1, int_type, &(Parallelization.ProcPEPsi), 1, critical);

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

  if(!ParseForParameter(verbose,FileBuffer,"DoOutOrbitals", 0, 1, 1, int_type, &(Switches.DoOutOrbitals), 1, optional)) {
    Switches.DoOutOrbitals = 0;
  } else {
    if (Switches.DoOutOrbitals < 0) Switches.DoOutOrbitals = 0;
    if (Switches.DoOutOrbitals > 1) Switches.DoOutOrbitals = 1;
  }
  ParseForParameter(verbose,FileBuffer,"DoOutVis", 0, 1, 1, int_type, &(Switches.DoOutVis), 1, critical);
  if (Switches.DoOutVis < 0) Switches.DoOutVis = 0;
  if (Switches.DoOutVis > 1) Switches.DoOutVis = 1;
  if (!ParseForParameter(verbose,FileBuffer,"VectorPlane", 0, 1, 1, int_type, &(LocalizedOrbitals.VectorPlane), 1, optional))
    LocalizedOrbitals.VectorPlane = -1;
  if (!ParseForParameter(verbose,FileBuffer,"VectorCut", 0, 1, 1, double_type, &(LocalizedOrbitals.VectorCut), 1, optional))
    LocalizedOrbitals.VectorCut = 0.;
  ParseForParameter(verbose,FileBuffer,"DoOutMes", 0, 1, 1, int_type, &(Switches.DoOutMes), 1, critical);
  if (Switches.DoOutMes < 0) Switches.DoOutMes = 0;
  if (Switches.DoOutMes > 1) Switches.DoOutMes = 1;
  if (!ParseForParameter(verbose,FileBuffer,"DoOutCurr", 0, 1, 1, int_type, &(Switches.DoOutCurrent), 1, optional))
    Switches.DoOutCurrent = 0;
  if (Switches.DoOutCurrent < 0) Switches.DoOutCurrent = 0;
  if (Switches.DoOutCurrent > 1) Switches.DoOutCurrent = 1;
  ParseForParameter(verbose,FileBuffer,"AddGramSch", 0, 1, 1, int_type, &(LocalizedOrbitals.UseAddGramSch), 1, critical);
  if (LocalizedOrbitals.UseAddGramSch < 0) LocalizedOrbitals.UseAddGramSch = 0;
  if (LocalizedOrbitals.UseAddGramSch > 2) LocalizedOrbitals.UseAddGramSch = 2;
  if(!ParseForParameter(verbose,FileBuffer,"DoWannier", 0, 1, 1, int_type, &(Switches.DoWannier), 1, optional)) {
    Switches.DoWannier = 0;
  } else {
    if (Switches.DoWannier < 0) Switches.DoWannier = 0;
    if (Switches.DoWannier > 1) Switches.DoWannier = 1;
  }
  if(!ParseForParameter(verbose,FileBuffer,"CommonWannier", 0, 1, 1, int_type, &(LocalizedOrbitals.CommonWannier), 1, optional)) {
    LocalizedOrbitals.CommonWannier = 0;
  } else {
    if (LocalizedOrbitals.CommonWannier < 0) LocalizedOrbitals.CommonWannier = 0;
    if (LocalizedOrbitals.CommonWannier > 4) LocalizedOrbitals.CommonWannier = 4;
  }
  if(!ParseForParameter(verbose,FileBuffer,"SawtoothStart", 0, 1, 1, double_type, &(LocalizedOrbitals.SawtoothStart), 1, optional)) {
    LocalizedOrbitals.SawtoothStart = 0.01;
  } else {
    if (LocalizedOrbitals.SawtoothStart < 0.) LocalizedOrbitals.SawtoothStart = 0.;
    if (LocalizedOrbitals.SawtoothStart > 1.) LocalizedOrbitals.SawtoothStart = 1.;
  }

  if (ParseForParameter(verbose,FileBuffer,"DoConstrainedMD", 0, 1, 1, int_type, &(Switches.DoConstrainedMD), 1, optional))
    if (Switches.DoConstrainedMD < 0)
      Switches.DoConstrainedMD = 0;
  ParseForParameter(verbose,FileBuffer,"MaxOuterStep", 0, 1, 1, int_type, &(StepCounts.MaxOuterStep), 1, critical);
  if (!ParseForParameter(verbose,FileBuffer,"Deltat", 0, 1, 1, double_type, &(Deltat), 1, optional))
    Deltat = 1;
  ParseForParameter(verbose,FileBuffer,"OutVisStep", 0, 1, 1, int_type, &(StepCounts.OutVisStep), 1, optional);
  ParseForParameter(verbose,FileBuffer,"OutSrcStep", 0, 1, 1, int_type, &(StepCounts.OutSrcStep), 1, optional);
  ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(World::getInstance().getThermostats()->TargetTemp), 1, optional);
  //ParseForParameter(verbose,FileBuffer,"Thermostat", 0, 1, 1, int_type, &(ScaleTempStep), 1, optional);
  if (!ParseForParameter(verbose,FileBuffer,"EpsWannier", 0, 1, 1, double_type, &(LocalizedOrbitals.EpsWannier), 1, optional))
    LocalizedOrbitals.EpsWannier = 1e-8;

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

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

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

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

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

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

  // IonsInitRead

  ParseForParameter(verbose,FileBuffer,"RCut", 0, 1, 1, double_type, &(PlaneWaveSpecifics.RCut), 1, critical);
  ParseForParameter(verbose,FileBuffer,"IsAngstroem", 0, 1, 1, int_type, &(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, &(RelativeCoord) , 1, optional))
    RelativeCoord = 0;
  if (!ParseForParameter(verbose,FileBuffer,"StructOpt", 0, 1, 1, int_type, &(StructOpt), 1, optional))
    StructOpt = 0;

  // 3. parse the molecule in
  molecule *mol = World::getInstance().createMolecule();
  LoadMolecule(mol, FileBuffer, World::getInstance().getPeriode(), FastParsing);
  //mol->SetNameFromFilename(filename);
  mol->ActiveFlag = true;
  //MolList->insert(mol);

  // 4. dissect the molecule into connected subgraphs
  // don't do this here ...
  //MolList->DissectMoleculeIntoConnectedSubgraphs(mol,this);
  //delete(mol);

  delete(FileBuffer);
}

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

/** Prints MaxTypes and list of elements to strea,
 * \param *file output stream
 * \param &allatoms vector of all atoms in the system, such as by World::getAllAtoms()
 * \param &ZtoIndexMap map of which atoms belong to which ion number
 */
void PcpParser::OutputElements(ostream *file, vector<atom *> &allatoms, map<int, int> &ZtoIndexMap)
{
  map<int, int> PresentElements;
  pair <   map<int, int>::iterator, bool > Inserter;
  // insert all found elements into the map
  for (vector<atom *>::iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
    Inserter = PresentElements.insert(pair<int, int>((*AtomRunner)->type->Z, 1));
    if (!Inserter.second) // increase if present
      Inserter.first->second += 1;
  }
  // print total element count
  *file << "MaxTypes\t" << PresentElements.size() <<  "\t# maximum number of different ion types" << endl;
  // print element list
  *file << "# Ion type data (PP = PseudoPotential, Z = atomic number)" << endl;
  *file << "#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\t# chemical name, symbol" << endl;
  // elements are due to map sorted by Z value automatically, hence just count through them
  int counter = 1;
  for(map<int, int>::const_iterator iter=PresentElements.begin(); iter!=PresentElements.end();++iter) {
    const element * const elemental = World::getInstance().getPeriode()->FindElement(iter->first);
    ZtoIndexMap.insert( pair<int,int> (iter->first, counter) );
    *file << "Ion_Type" << counter++ << "\t" << iter->second << "\t" << elemental->Z << "\t1.0\t3\t3\t" << fixed << setprecision(11) << showpoint << elemental->mass << "\t" << elemental->name << "\t" << elemental->symbol <<endl;
  }
}

/** Output all atoms one per line.
 * \param *file output stream
 * \param &allatoms vector of all atoms in the system, such as by World::getAllAtoms()
 * \param &ZtoIndexMap map of which atoms belong to which ion number
 */
void PcpParser::OutputAtoms(ostream *file, vector<atom *> &allatoms, map<int, int> &ZtoIndexMap)
{
  *file << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
  map<int, int> ZtoCountMap;
  pair <   map<int, int>::iterator, bool > Inserter;
  for (vector<atom *>::iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
    Inserter = ZtoCountMap.insert( pair<int, int>((*AtomRunner)->type->Z, 1) );
    if (!Inserter.second)
      Inserter.first->second += 1;
    (*AtomRunner)->OutputArrayIndexed(file, &ZtoIndexMap[(*AtomRunner)->type->Z], &ZtoCountMap[Inserter.first->second], NULL);
  }
}

/** Reading of Thermostat related values from parameter file.
 * \param *fb file buffer containing the config file
 */
void PcpParser::ParseThermostats(class ConfigFileBuffer * const fb)
{
  char * const thermo = new char[12];
  const int verbose = 0;
  class 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)) {
    if (strcmp(thermo, Thermostats->ThermostatNames[0]) == 0) { // None
      if (Thermostats->ThermostatImplemented[0] == 1) {
        Thermostats->Thermostat = None;
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
        Thermostats->Thermostat = None;
      }
    } else if (strcmp(thermo, Thermostats->ThermostatNames[1]) == 0) { // Woodcock
      if (Thermostats->ThermostatImplemented[1] == 1) {
        Thermostats->Thermostat = Woodcock;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read scaling frequency
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
        Thermostats->Thermostat = None;
      }
    } else if (strcmp(thermo, Thermostats->ThermostatNames[2]) == 0) { // Gaussian
      if (Thermostats->ThermostatImplemented[2] == 1) {
        Thermostats->Thermostat = Gaussian;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read collision rate
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
        Thermostats->Thermostat = None;
      }
    } else if (strcmp(thermo, Thermostats->ThermostatNames[3]) == 0) { // Langevin
      if (Thermostats->ThermostatImplemented[3] == 1) {
        Thermostats->Thermostat = Langevin;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read gamma
        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &Thermostats->alpha, 1, optional)) {
          DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << Thermostats->alpha << "." << endl);
        } else {
          Thermostats->alpha = 1.;
        }
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
        Thermostats->Thermostat = None;
      }
    } else if (strcmp(thermo, Thermostats->ThermostatNames[4]) == 0) { // Berendsen
      if (Thermostats->ThermostatImplemented[4] == 1) {
        Thermostats->Thermostat = Berendsen;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read \tau_T
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
        Thermostats->Thermostat = None;
      }
    } else if (strcmp(thermo, Thermostats->ThermostatNames[5]) == 0) { // Nose-Hoover
      if (Thermostats->ThermostatImplemented[5] == 1) {
        Thermostats->Thermostat = NoseHoover;
        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->HooverMass, 1, critical); // read Hoovermass
        Thermostats->alpha = 0.;
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
        Thermostats->Thermostat = None;
      }
    } else {
      DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
      Thermostats->Thermostat = None;
    }
  } else {
    if ((Thermostats->TargetTemp != 0))
      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    Thermostats->Thermostat = None;
  }
  delete[](thermo);
};