/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include <iostream>
#include <iomanip>

#include "Atom/atom.hpp"
#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "PcpParser.hpp"
#include "Parser/ConfigFileBuffer.hpp"
#include "Parser/PcpParser_helper.hpp"
#include "Thermostats/ThermoStatContainer.hpp"
#include "World.hpp"


// declare specialized static variables
const std::string FormatParserTrait<pcp>::name = "pcp";
const std::string FormatParserTrait<pcp>::suffix = "conf";
const ParserTypes FormatParserTrait<pcp>::type = pcp;

FormatParser< pcp >::StructParallelization::StructParallelization() :
  ProcPEGamma(8),
  ProcPEPsi(1)
{}

FormatParser< pcp >::StructParallelization::~StructParallelization()
{}

FormatParser< pcp >::StructPaths::StructPaths() :
  databasepath(NULL),
  configname(NULL),
  mainname(NULL),
  defaultpath(NULL),
  pseudopotpath(NULL)
{}

FormatParser< pcp >::StructPaths::~StructPaths()
{}

FormatParser< pcp >::StructSwitches::StructSwitches() :
  DoConstrainedMD(0),
  DoOutVis(0),
  DoOutMes(1),
  DoOutNICS(0),
  DoOutOrbitals(0),
  DoOutCurrent(0),
  DoFullCurrent(0),
  DoPerturbation(0),
  DoWannier(0)
{}

FormatParser< pcp >::StructSwitches::~StructSwitches()
{}

FormatParser< pcp >::StructLocalizedOrbitals::StructLocalizedOrbitals() :
  CommonWannier(0),
  SawtoothStart(0.01),
  VectorPlane(0),
  VectorCut(0),
  UseAddGramSch(1),
  Seed(1),
  EpsWannier(1e-7)
{}

FormatParser< pcp >::StructLocalizedOrbitals::~StructLocalizedOrbitals()
{}

FormatParser< pcp >::StructStepCounts::StructStepCounts() :
  MaxMinStopStep(1),
  InitMaxMinStopStep(1),
  OutVisStep(10),
  OutSrcStep(5),
  MaxPsiStep(0),
  MaxOuterStep(0),
  MaxMinStep(100),
  RelEpsTotalEnergy(1e-07),
  RelEpsKineticEnergy(1e-05),
  MaxMinGapStopStep(0),
  MaxInitMinStep(100),
  InitRelEpsTotalEnergy(1e-05),
  InitRelEpsKineticEnergy(0.0001),
  InitMaxMinGapStopStep(0)
{}

FormatParser< pcp >::StructStepCounts::~StructStepCounts()
{}

FormatParser< pcp >::StructPlaneWaveSpecifics::StructPlaneWaveSpecifics() :
  PsiType(0),
  MaxPsiDouble(0),
  PsiMaxNoUp(0),
  PsiMaxNoDown(0),
  ECut(128),
  MaxLevel(5),
  RiemannTensor(0),
  LevRFactor(0),
  RiemannLevel(0),
  Lev0Factor(2),
  RTActualUse(0),
  AddPsis(0),
  RCut(20)
{}

FormatParser< pcp >::StructPlaneWaveSpecifics::~StructPlaneWaveSpecifics()
{}

/** Constructor of PcpParser.
 *
 */
FormatParser< pcp >::FormatParser() :
  FormatParser_common(NULL),
  FastParsing(false),
  Deltat(0.01),
  IsAngstroem(1),
  RelativeCoord(0),
  StructOpt(0),
  MaxTypes(0)
{}

/** Destructor of PcpParser.
 *
 */
FormatParser< pcp >::~FormatParser()
{}

void FormatParser< pcp >::load(std::istream* file)
{
  if (file->fail()) {
    ELOG(1, "could not access given file");
    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
  if (!ParseForParameter(verbose,FileBuffer, "mainname", 0, 1, 1, string_type, (Paths.mainname), 1, critical)) {
    ELOG(1, "mainname is missing, is file empty?");
  } else {
    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 *cell_size = new double[6];
    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];
    World::getInstance().setDomain(cell_size);
    delete[] cell_size;
    //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,"PsiMaxNoUp", 0, 1, 1, int_type, &(PlaneWaveSpecifics.PsiMaxNoUp), 1, optional);
      ParseForParameter(verbose,FileBuffer,"PsiMaxNoDown", 0, 1, 1, int_type, &(PlaneWaveSpecifics.PsiMaxNoDown), 1, optional);
      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,"MaxPsiDouble", 0, 1, 1, int_type, &(PlaneWaveSpecifics.MaxPsiDouble), 1, optional);
      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();
  mol->ActiveFlag = true;
  // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
  World::getInstance().getMolecules()->insert(mol);
  LoadMolecule(mol, FileBuffer, World::getInstance().getPeriode(), FastParsing);

  // refresh atom::nr and atom::name
  mol->getAtomCount();

  // 4. dissect the molecule into connected subgraphs
  // don't do this here ...
  //FragmentationSubgraphDissection();
  //delete(mol);

  delete(FileBuffer);
}

/**
 * Saves the \a atoms into as a PCP file.
 *
 * \param file where to save the state
 * \param atoms atoms to store
 */
void FormatParser< pcp >::save(std::ostream* file, const std::vector<atom *> &atoms)
{
  LOG(0, "Saving changes to pcp.");

  const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
  ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
  if (!file->fail()) {
    // calculate number of Psis
    CalculateOrbitals(atoms);
    *file << "# ParallelCarParinello - main configuration file - created with molecuilder" << endl;
    *file << endl;
    if (Paths.mainname != NULL)
      *file << "mainname\t" << Paths.mainname << "\t# programm name (for runtime files)" << endl;
    else
      *file << "mainname\tpcp\t# programm name (for runtime files)" << endl;
    if (Paths.defaultpath != NULL)
      *file << "defaultpath\t" << Paths.defaultpath << "\t# where to put files during runtime" << endl;
    else
      *file << "defaultpath\tnot specified\t# where to put files during runtime" << endl;
    if (Paths.pseudopotpath != NULL)
      *file << "pseudopotpath\t" << Paths.pseudopotpath << "\t# where to find pseudopotentials" << endl;
    else
      *file << "pseudopotpath\tnot specified\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;
    ASSERT(Thermostats != NULL, "FormatParser< pcp >::save() - Thermostats not initialized!");
    *file << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
    *file << Thermostats->activeThermostat->writeParams();
    *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 << domain.at(0,0) << "\t" << endl;
    *file << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
    *file << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\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;
    *file << "MaxPsiDouble\t" << PlaneWaveSpecifics.MaxPsiDouble << "\t# here: specifying both maximum number of SpinUp- and -Down-states" << endl;
    *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;
    *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;
    map<int, int> ZtoIndexMap;
    OutputElements(file, atoms, ZtoIndexMap);
    OutputAtoms(file, atoms, ZtoIndexMap);
  } else {
    ELOG(1, "Cannot open output file.");
  }
}


/** Counts necessary number of valence electrons and returns number and SpinType.
 * \param &allatoms all atoms to store away
 */
void FormatParser< pcp >::CalculateOrbitals(const std::vector<atom *> &allatoms)
{
  PlaneWaveSpecifics.MaxPsiDouble = PlaneWaveSpecifics.PsiMaxNoDown = PlaneWaveSpecifics.PsiMaxNoUp = PlaneWaveSpecifics.PsiType = 0;
  for (vector<atom *>::const_iterator runner = allatoms.begin(); runner != allatoms.end(); ++runner) {
    PlaneWaveSpecifics.MaxPsiDouble += (*runner)->getType()->getNoValenceOrbitals();
  }
  cout << PlaneWaveSpecifics.MaxPsiDouble << endl;
  PlaneWaveSpecifics.PsiMaxNoDown = PlaneWaveSpecifics.MaxPsiDouble/2 + (PlaneWaveSpecifics.MaxPsiDouble % 2);
  PlaneWaveSpecifics.PsiMaxNoUp = PlaneWaveSpecifics.MaxPsiDouble/2;
  PlaneWaveSpecifics.MaxPsiDouble /= 2;
  PlaneWaveSpecifics.PsiType = (PlaneWaveSpecifics.PsiMaxNoDown == PlaneWaveSpecifics.PsiMaxNoUp) ? 0 : 1;
  if ((PlaneWaveSpecifics.PsiType == 1) && (Parallelization.ProcPEPsi < 2) && ((PlaneWaveSpecifics.PsiMaxNoDown != 1) || (PlaneWaveSpecifics.PsiMaxNoUp != 0))) {
    Parallelization.ProcPEGamma /= 2;
    Parallelization.ProcPEPsi *= 2;
  } else {
    Parallelization.ProcPEGamma *= Parallelization.ProcPEPsi;
    Parallelization.ProcPEPsi = 1;
  }
  cout << PlaneWaveSpecifics.PsiMaxNoDown << ">" << PlaneWaveSpecifics.PsiMaxNoUp << endl;
  if (PlaneWaveSpecifics.PsiMaxNoDown > PlaneWaveSpecifics.PsiMaxNoUp) {
    StepCounts.InitMaxMinStopStep = StepCounts.MaxMinStopStep = PlaneWaveSpecifics.PsiMaxNoDown;
    cout << PlaneWaveSpecifics.PsiMaxNoDown << " " << StepCounts.InitMaxMinStopStep << endl;
  } else {
    StepCounts.InitMaxMinStopStep = StepCounts.MaxMinStopStep = PlaneWaveSpecifics.PsiMaxNoUp;
    cout << PlaneWaveSpecifics.PsiMaxNoUp << " " << StepCounts.InitMaxMinStopStep << 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 FormatParser< pcp >::OutputElements(ostream *file, const std::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 *>::const_iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
    Inserter = PresentElements.insert(pair<int, int>((*AtomRunner)->getType()->getAtomicNumber(), 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;
  *file << 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->getAtomicNumber() << "\t1.0\t3\t3\t" << fixed << setprecision(11) << showpoint << elemental->getMass() << "\t" << elemental->getName() << "\t" << elemental->getSymbol() <<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 FormatParser< pcp >::OutputAtoms(ostream *file, const std::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;
  map<atom *, int> AtomtoCountMap;
  pair <   map<int, int>::iterator, bool > Inserter;
  bool ContinueStatus = true;
  bool AddNewLine = false;
  size_t step = 0;
  do {
    int nr = 0;
    ContinueStatus = false;
    for (vector<atom *>::const_iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
      if ((*AtomRunner)->getTrajectorySize() > step) {
        if (step == 0) {  // fill list only on first step
          Inserter = ZtoCountMap.insert( pair<int, int>((*AtomRunner)->getType()->getAtomicNumber(), 1) );
          if (!Inserter.second)
            Inserter.first->second += 1;
          AtomtoCountMap.insert( make_pair((*AtomRunner), Inserter.first->second) );
        }
        if (AddNewLine) {
          AddNewLine = false;
          *file << endl;
        }
        const int Z = (*AtomRunner)->getType()->getAtomicNumber();
        *file << "Ion_Type" << ZtoIndexMap[Z] << "_" << AtomtoCountMap[(*AtomRunner)] << "\t"  << fixed << setprecision(9) << showpoint;
        *file << (*AtomRunner)->atStep(0, step)
            << "\t" << (*AtomRunner)->atStep(1,step)
            << "\t" << (*AtomRunner)->atStep(2,step);
        *file << "\t" << (int)((*AtomRunner)->getFixedIon());
        if ((*AtomRunner)->getAtomicVelocityAtStep(step).Norm() > MYEPSILON)
          *file << "\t" << scientific << setprecision(6)
              << (*AtomRunner)->getAtomicVelocityAtStep(step)[0] << "\t"
              << (*AtomRunner)->getAtomicVelocityAtStep(step)[1] << "\t"
              << (*AtomRunner)->getAtomicVelocityAtStep(step)[2] << "\t";
        *file << " # molecule nr " << nr++ << endl;
        ContinueStatus = true;  // as long as a single atom still has trajectory points, continue
      }
    }
    ++step;
    if (ContinueStatus)
      AddNewLine = true;
  } while (ContinueStatus);
}

/** Reading of Thermostat related values from parameter file.
 * \param *fb file buffer containing the config file
 */
void FormatParser< pcp >::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)) {
    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);
};

bool FormatParser< pcp >::operator==(const FormatParser< pcp >& b) const
{
  ASSERT(Parallelization.ProcPEGamma == b.Parallelization.ProcPEGamma, "PcpParser ==: ProcPEGamma not");
  ASSERT(Parallelization.ProcPEPsi == b.Parallelization.ProcPEPsi, "PcpParser ==: ProcPEPsi not");

  if ((Paths.databasepath != NULL) && (b.Paths.databasepath != NULL))
    ASSERT(strcmp(Paths.databasepath, b.Paths.databasepath), "PcpParser ==: databasepath not");
  if ((Paths.configname != NULL) && (b.Paths.configname != NULL))
    ASSERT(strcmp(Paths.configname, b.Paths.configname), "PcpParser ==: configname not");
  if ((Paths.mainname != NULL) && (b.Paths.mainname != NULL))
    ASSERT(strcmp(Paths.mainname, b.Paths.mainname), "PcpParser ==: mainname not");
  if ((Paths.defaultpath != NULL) && (b.Paths.defaultpath != NULL))
    ASSERT(strcmp(Paths.defaultpath, b.Paths.defaultpath), "PcpParser ==: defaultpath not");
  if ((Paths.pseudopotpath != NULL) && (b.Paths.pseudopotpath != NULL))
    ASSERT(strcmp(Paths.pseudopotpath, b.Paths.pseudopotpath), "PcpParser ==: pseudopotpath not");

  ASSERT(Switches.DoConstrainedMD == b.Switches.DoConstrainedMD, "PcpParser ==: DoConstrainedMD not");
  ASSERT(Switches.DoOutVis == b.Switches.DoOutVis, "PcpParser ==: DoOutVis not");
  ASSERT(Switches.DoOutMes == b.Switches.DoOutMes, "PcpParser ==: DoOutMes not");
  ASSERT(Switches.DoOutNICS == b.Switches.DoOutNICS, "PcpParser ==: DoOutNICS not");
  ASSERT(Switches.DoOutOrbitals == b.Switches.DoOutOrbitals, "PcpParser ==: DoOutOrbitals not");
  ASSERT(Switches.DoOutCurrent == b.Switches.DoOutCurrent, "PcpParser ==: DoOutCurrent not");
  ASSERT(Switches.DoFullCurrent == b.Switches.DoFullCurrent, "PcpParser ==: DoFullCurrent not");
  ASSERT(Switches.DoPerturbation == b.Switches.DoPerturbation, "PcpParser ==: DoPerturbation not");
  ASSERT(Switches.DoWannier == b.Switches.DoWannier, "PcpParser ==: DoWannier not");

  ASSERT(LocalizedOrbitals.CommonWannier == b.LocalizedOrbitals.CommonWannier, "PcpParser ==: CommonWannier not");
  ASSERT(LocalizedOrbitals.SawtoothStart == b.LocalizedOrbitals.SawtoothStart, "PcpParser ==: SawtoothStart not");
  ASSERT(LocalizedOrbitals.VectorPlane == b.LocalizedOrbitals.VectorPlane, "PcpParser ==: VectorPlane not");
  ASSERT(LocalizedOrbitals.VectorCut == b.LocalizedOrbitals.VectorCut, "PcpParser ==: VectorCut not");
  ASSERT(LocalizedOrbitals.UseAddGramSch == b.LocalizedOrbitals.UseAddGramSch, "PcpParser ==: UseAddGramSch not");
  ASSERT(LocalizedOrbitals.Seed == b.LocalizedOrbitals.Seed, "PcpParser ==: Seed not");
  ASSERT(LocalizedOrbitals.EpsWannier == b.LocalizedOrbitals.EpsWannier, "PcpParser ==: EpsWannier not");

  ASSERT(StepCounts.MaxMinStopStep == b.StepCounts.MaxMinStopStep, "PcpParser ==: MaxMinStopStep not");
  ASSERT(StepCounts.InitMaxMinStopStep == b.StepCounts.InitMaxMinStopStep, "PcpParser ==: InitMaxMinStopStep not");
  ASSERT(StepCounts.OutVisStep == b.StepCounts.OutVisStep, "PcpParser ==: OutVisStep not");
  ASSERT(StepCounts.OutSrcStep == b.StepCounts.OutSrcStep, "PcpParser ==: OutSrcStep not");
  ASSERT(StepCounts.MaxPsiStep == b.StepCounts.MaxPsiStep, "PcpParser ==: MaxPsiStep not");
  ASSERT(StepCounts.MaxOuterStep == b.StepCounts.MaxOuterStep, "PcpParser ==: MaxOuterStep not");
  ASSERT(StepCounts.MaxMinStep == b.StepCounts.MaxMinStep, "PcpParser ==: MaxMinStep not");
  ASSERT(StepCounts.RelEpsTotalEnergy == b.StepCounts.RelEpsTotalEnergy, "PcpParser ==: RelEpsTotalEnergy not");
  ASSERT(StepCounts.MaxMinGapStopStep == b.StepCounts.MaxMinGapStopStep, "PcpParser ==: MaxMinGapStopStep not");
  ASSERT(StepCounts.MaxInitMinStep == b.StepCounts.MaxInitMinStep, "PcpParser ==: MaxInitMinStep not");
  ASSERT(StepCounts.InitRelEpsTotalEnergy == b.StepCounts.InitRelEpsTotalEnergy, "PcpParser ==: InitRelEpsTotalEnergy not");
  ASSERT(StepCounts.InitRelEpsKineticEnergy == b.StepCounts.InitRelEpsKineticEnergy, "PcpParser ==: InitRelEpsKineticEnergy not");
  ASSERT(StepCounts.InitMaxMinGapStopStep == b.StepCounts.InitMaxMinGapStopStep, "PcpParser ==: InitMaxMinGapStopStep not");

  ASSERT(PlaneWaveSpecifics.PsiType == b.PlaneWaveSpecifics.PsiType, "PcpParser ==: PsiType not");
  ASSERT(PlaneWaveSpecifics.MaxPsiDouble == b.PlaneWaveSpecifics.MaxPsiDouble, "PcpParser ==: MaxPsiDouble not");
  ASSERT(PlaneWaveSpecifics.PsiMaxNoUp == b.PlaneWaveSpecifics.PsiMaxNoUp, "PcpParser ==: PsiMaxNoUp not");
  ASSERT(PlaneWaveSpecifics.PsiMaxNoDown == b.PlaneWaveSpecifics.PsiMaxNoDown, "PcpParser ==: PsiMaxNoDown not");
  ASSERT(PlaneWaveSpecifics.ECut == b.PlaneWaveSpecifics.ECut, "PcpParser ==: ECut not");
  ASSERT(PlaneWaveSpecifics.MaxLevel == b.PlaneWaveSpecifics.MaxLevel, "PcpParser ==: MaxLevel not");
  ASSERT(PlaneWaveSpecifics.RiemannTensor == b.PlaneWaveSpecifics.RiemannTensor, "PcpParser ==: RiemannTensor not");
  ASSERT(PlaneWaveSpecifics.LevRFactor == b.PlaneWaveSpecifics.LevRFactor, "PcpParser ==: LevRFactor not");
  ASSERT(PlaneWaveSpecifics.RiemannLevel == b.PlaneWaveSpecifics.RiemannLevel, "PcpParser ==: RiemannLevel not");
  ASSERT(PlaneWaveSpecifics.Lev0Factor == b.PlaneWaveSpecifics.Lev0Factor, "PcpParser ==: Lev0Factor not");
  ASSERT(PlaneWaveSpecifics.RTActualUse == b.PlaneWaveSpecifics.RTActualUse, "PcpParser ==: RTActualUse not");
  ASSERT(PlaneWaveSpecifics.AddPsis == b.PlaneWaveSpecifics.AddPsis, "PcpParser ==: AddPsis not");
  ASSERT(PlaneWaveSpecifics.AddPsis == b.PlaneWaveSpecifics.AddPsis, "PcpParser ==: AddPsis not");
  ASSERT(PlaneWaveSpecifics.RCut == b.PlaneWaveSpecifics.RCut, "PcpParser ==: RCut not");

  ASSERT(FastParsing == b.FastParsing, "PcpParser ==: FastParsing not");

  ASSERT(Deltat == b.Deltat, "PcpParser ==: Deltat not");
  ASSERT(IsAngstroem == b.IsAngstroem, "PcpParser ==: IsAngstroem not");
  ASSERT(RelativeCoord == b.RelativeCoord, "PcpParser ==: RelativeCoord not");
  ASSERT(StructOpt == b.StructOpt, "PcpParser ==: StructOpt not");
  ASSERT(MaxTypes == b.MaxTypes, "PcpParser ==: MaxTypes not");
  ASSERT(basis == b.basis, "PcpParser ==: basis not");

  return true;
}