source: src/Parser/PcpParser.cpp@ 1a6bda

/*
 * 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.
 */

/*
 * 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 "Actions/FragmentationAction/SubgraphDissectionAction.hpp"
#include "atom.hpp"
#include "config.hpp"
#include "ConfigFileBuffer.hpp"
#include "element.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "molecule.hpp"
#include "PcpParser.hpp"
#include "periodentafel.hpp"
#include "ThermoStatContainer.hpp"
#include "World.hpp"
#include "Box.hpp"


PcpParser::StructParallelization::StructParallelization() :
  ProcPEGamma(8),
  ProcPEPsi(1)
{}

PcpParser::StructParallelization::~StructParallelization()
{}

PcpParser::StructPaths::StructPaths() :
  databasepath(NULL),
  configname(NULL),
  mainname(NULL),
  defaultpath(NULL),
  pseudopotpath(NULL)
{}

PcpParser::StructPaths::~StructPaths()
{}

PcpParser::StructSwitches::StructSwitches() :
  DoConstrainedMD(0),
  DoOutVis(0),
  DoOutMes(1),
  DoOutNICS(0),
  DoOutOrbitals(0),
  DoOutCurrent(0),
  DoFullCurrent(0),
  DoPerturbation(0),
  DoWannier(0)
{}

PcpParser::StructSwitches::~StructSwitches()
{}

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

PcpParser::StructLocalizedOrbitals::~StructLocalizedOrbitals()
{}

PcpParser::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)
{}

PcpParser::StructStepCounts::~StructStepCounts()
{}

PcpParser::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)
{}

PcpParser::StructPlaneWaveSpecifics::~StructPlaneWaveSpecifics()
{}

/** Constructor of PcpParser.
 *
 */
PcpParser::PcpParser() :
  FastParsing(false),
  Deltat(0.01),
  IsAngstroem(1),
  RelativeCoord(0),
  StructOpt(0),
  MaxTypes(0)
{}

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

  // 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 PcpParser::save(std::ostream* file, const std::vector<atom *> &atoms)
{
  DoLog(0) && (Log() << Verbose(0) << "Saving changes to pcp." << std::endl);

  const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
  class 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, "PcpParser::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 {
    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open output file." << endl);
  }
}


/** Counts necessary number of valence electrons and returns number and SpinType.
 * \param &allatoms all atoms to store away
 */
void PcpParser::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 PcpParser::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 PcpParser::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;
  pair <   map<int, int>::iterator, bool > Inserter;
  int nr = 0;
  for (vector<atom *>::const_iterator AtomRunner = allatoms.begin();AtomRunner != allatoms.end();++AtomRunner) {
    Inserter = ZtoCountMap.insert( pair<int, int>((*AtomRunner)->getType()->getAtomicNumber(), 1) );
    if (!Inserter.second)
      Inserter.first->second += 1;
    const int Z = (*AtomRunner)->getType()->getAtomicNumber();
    *file << "Ion_Type" << ZtoIndexMap[Z] << "_" << ZtoCountMap[Z] << "\t"  << fixed << setprecision(9) << showpoint;
    *file << (*AtomRunner)->at(0) << "\t" << (*AtomRunner)->at(1) << "\t" << (*AtomRunner)->at(2);
    *file << "\t" << (*AtomRunner)->FixedIon;
    if ((*AtomRunner)->AtomicVelocity.Norm() > MYEPSILON)
      *file << "\t" << scientific << setprecision(6) << (*AtomRunner)->AtomicVelocity[0] << "\t" << (*AtomRunner)->AtomicVelocity[1] << "\t" << (*AtomRunner)->AtomicVelocity[2] << "\t";
    *file << " # molecule nr " << nr++ << endl;
  }
}

/** 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)) {
    Thermostats->makeActive(thermo,fb);
  } else {
    if ((Thermostats->TargetTemp != 0))
      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    Thermostats->chooseNone();
  }
  delete[](thermo);
};

bool PcpParser::operator==(const PcpParser& 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;
}