Index: src/Actions/CmdAction/BondLengthTableAction.cpp
===================================================================
--- src/Actions/CmdAction/BondLengthTableAction.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Actions/CmdAction/BondLengthTableAction.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -9,4 +9,5 @@
 
 #include "Actions/CmdAction/BondLengthTableAction.hpp"
+#include "bondgraph.hpp"
 #include "config.hpp"
 #include "log.hpp"
Index: src/Actions/FragmentationAction/DepthFirstSearchAction.cpp
===================================================================
--- src/Actions/FragmentationAction/DepthFirstSearchAction.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Actions/FragmentationAction/DepthFirstSearchAction.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -10,4 +10,5 @@
 #include "Actions/FragmentationAction/DepthFirstSearchAction.hpp"
 #include "atom.hpp"
+#include "bondgraph.hpp"
 #include "config.hpp"
 #include "log.hpp"
Index: src/Actions/FragmentationAction/FragmentationAction.cpp
===================================================================
--- src/Actions/FragmentationAction/FragmentationAction.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Actions/FragmentationAction/FragmentationAction.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -10,4 +10,5 @@
 #include "Actions/FragmentationAction/FragmentationAction.hpp"
 #include "atom.hpp"
+#include "bondgraph.hpp"
 #include "config.hpp"
 #include "log.hpp"
Index: src/Actions/MapOfActions.cpp
===================================================================
--- src/Actions/MapOfActions.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Actions/MapOfActions.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -266,5 +266,5 @@
   generic.insert("nonconvex-envelope");
 	generic.insert("pair-correlation");
-//	generic.insert("parse-xyz");
+	generic.insert("parse-xyz");
 //  generic.insert("principal-axis-system");
   generic.insert("remove-atom");
Index: src/Actions/ParserAction/LoadXyzAction.cpp
===================================================================
--- src/Actions/ParserAction/LoadXyzAction.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Actions/ParserAction/LoadXyzAction.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -8,11 +8,14 @@
 #include "Helpers/MemDebug.hpp"
 
+using namespace std;
+
 #include "Actions/ParserAction/LoadXyzAction.hpp"
 #include "Parser/XyzParser.hpp"
 
 #include <iostream>
+#include <set>
 #include <string>
+#include <vector>
 
-using namespace std;
 
 #include "UIElements/UIFactory.hpp"
@@ -21,5 +24,8 @@
 
 #include "atom.hpp"
+#include "log.hpp"
 #include "molecule.hpp"
+#include "verbose.hpp"
+#include "World.hpp"
 
 /****** ParserLoadXyzAction *****/
@@ -37,5 +43,5 @@
 //};
 
-const char ParserLoadXyzAction::NAME[] = "LoadXyz";
+const char ParserLoadXyzAction::NAME[] = "parse-xyz";
 
 ParserLoadXyzAction::ParserLoadXyzAction() :
@@ -48,15 +54,50 @@
 Action::state_ptr ParserLoadXyzAction::performCall() {
   string filename;
-  XyzParser parser;
   Dialog *dialog = UIFactory::getInstance().makeDialog();
 
-  dialog->queryString("filename",&filename, "Filename of the xyz file");
+  dialog->queryString(NAME,&filename, NAME);
 
   if(dialog->display()) {
+    DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl);
     // parse xyz file
     ifstream input;
     input.open(filename.c_str());
-    if (!input.fail())
+    if (!input.fail()) {
+      // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
+      set <atom*> UniqueList;
+      {
+        vector<atom *> ListBefore = World::getInstance().getAllAtoms();
+        for (vector<atom *>::iterator runner = ListBefore.begin();runner != ListBefore.end(); ++runner)
+          UniqueList.insert(*runner);
+      }
+      XyzParser parser; // briefly instantiate a parser which is removed at end of focus
       parser.load(&input);
+      {
+        vector<atom *> ListAfter = World::getInstance().getAllAtoms();
+        pair< set<atom *>::iterator, bool > Inserter;
+        if (UniqueList.size() != ListAfter.size()) { // only create if new atoms have been parsed
+          MoleculeListClass *molecules = World::getInstance().getMolecules();
+          molecule *mol= NULL;
+          if (molecules->ListOfMolecules.empty()) {
+            mol = World::getInstance().createMolecule();
+            molecules->insert(mol);
+          } else {
+            mol = *(molecules->ListOfMolecules.begin());
+          }
+          for (vector<atom *>::iterator runner = ListAfter.begin(); runner != ListAfter.end(); ++runner) {
+            Inserter = UniqueList.insert(*runner);
+            if (Inserter.second) { // if not present, then new (just parsed) atom, add ...
+              cout << "Adding new atom " << **runner << " to new mol." << endl;
+              mol->AddAtom(*runner);
+            }
+          }
+          mol->doCountAtoms();
+        } else {
+          cout << "No atoms parsed?" << endl;
+        }
+      }
+    } else {
+      DoeLog(1) && (eLog() << Verbose(1) << "Could not open file " << filename << "." << endl);
+    }
     input.close();
   }
Index: src/ConfigFileBuffer.cpp
===================================================================
--- src/ConfigFileBuffer.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/ConfigFileBuffer.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,182 @@
+/*
+ * ConfigFileBuffer.cpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#include "ConfigFileBuffer.hpp"
+#include "helpers.hpp"
+#include "lists.hpp"
+#include "log.hpp"
+#include "World.hpp"
+
+/******************************** Functions for class ConfigFileBuffer **********************/
+
+/** Structure containing compare function for Ion_Type sorting.
+ */
+struct IonTypeCompare {
+  bool operator()(const char* s1, const char *s2) const {
+    char number1[8];
+    char number2[8];
+    const char *dummy1, *dummy2;
+    //Log() << Verbose(0) << s1 << "  " << s2 << endl;
+    dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
+    dummy2 = strchr(dummy1, '_');
+    strncpy(number1, dummy1, dummy2-dummy1); // copy the number
+    number1[dummy2-dummy1]='\0';
+    dummy1 = strchr(s2, '_')+sizeof(char)*5;  // go just after "Ion_Type"
+    dummy2 = strchr(dummy1, '_');
+    strncpy(number2, dummy1, dummy2-dummy1); // copy the number
+    number2[dummy2-dummy1]='\0';
+    if (atoi(number1) != atoi(number2))
+      return (atoi(number1) < atoi(number2));
+    else {
+      dummy1 = strchr(s1, '_')+sizeof(char);
+      dummy1 = strchr(dummy1, '_')+sizeof(char);
+      dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
+      strncpy(number1, dummy1, dummy2-dummy1); // copy the number
+      number1[dummy2-dummy1]='\0';
+      dummy1 = strchr(s2, '_')+sizeof(char);
+      dummy1 = strchr(dummy1, '_')+sizeof(char);
+      dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
+      strncpy(number2, dummy1, dummy2-dummy1); // copy the number
+      number2[dummy2-dummy1]='\0';
+      return (atoi(number1) < atoi(number2));
+    }
+  }
+};
+
+/** Constructor for ConfigFileBuffer class.
+ */
+ConfigFileBuffer::ConfigFileBuffer() : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
+{
+};
+
+/** Constructor for ConfigFileBuffer class with filename to be parsed.
+ * \param *filename file name
+ */
+ConfigFileBuffer::ConfigFileBuffer(const char * const filename) : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
+{
+  InitFileBuffer(filename);
+}
+
+void ConfigFileBuffer::InitFileBuffer(const char * const filename)
+{
+  ifstream *file= new ifstream(filename);
+  InitFileBuffer(file);
+}
+
+void ConfigFileBuffer::InitFileBuffer(istream *file)
+{
+  char line[MAXSTRINGSIZE];
+
+  RemoveMapping();
+
+  // prescan number of lines
+  if (file->fail()) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "config file missing!" << endl);
+    return;
+  }
+  NoLines = 0; // we're overcounting by one
+  long file_position = file->tellg(); // mark current position
+  do {
+    file->getline(line, MAXSTRINGSIZE-1);
+    NoLines++;
+  } while (!file->eof());
+  file->clear();
+  file->seekg(file_position, ios::beg);
+  DoLog(1) && (Log() << Verbose(1) << NoLines-1 << " lines were recognized." << endl);
+
+  // allocate buffer's 1st dimension
+  if (buffer != NULL) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "FileBuffer->buffer is not NULL!" << endl);
+    return;
+  } else
+    buffer = new char *[NoLines];
+
+  // scan each line and put into buffer
+  int lines=0;
+  int i;
+  do {
+    buffer[lines] = new char[MAXSTRINGSIZE];
+    file->getline(buffer[lines], MAXSTRINGSIZE-1);
+    i = strlen(buffer[lines]);
+    buffer[lines][i] = '\n';
+    buffer[lines][i+1] = '\0';
+    lines++;
+  } while((!file->eof()) && (lines < NoLines));
+  DoLog(1) && (Log() << Verbose(1) << lines-1 << " lines were read into the buffer." << endl);
+  file->clear();
+  file->seekg(file_position, ios::beg);
+
+  InitMapping();
+}
+
+/** Destructor for ConfigFileBuffer class.
+ */
+ConfigFileBuffer::~ConfigFileBuffer()
+{
+  RemoveBuffer();
+  RemoveMapping();
+}
+
+
+/** Create trivial mapping.
+ */
+void ConfigFileBuffer::InitMapping()
+{
+  LineMapping = new int[NoLines];
+  for (int i=0;i<NoLines;i++)
+    LineMapping[i] = i;
+  MappingAllocated = true;
+}
+
+/** Remove allocated mapping.
+ */
+void ConfigFileBuffer::RemoveMapping()
+{
+  delete[](LineMapping);
+  MappingAllocated = false;
+}
+
+/** Remove allocated mapping.
+ */
+void ConfigFileBuffer::RemoveBuffer()
+{
+  for(int i=0;i<NoLines;++i)
+    delete[](buffer[i]);
+  delete[](buffer);
+}
+
+
+/** Creates a mapping for the \a *FileBuffer's lines containing the Ion_Type keyword such that they are sorted.
+ * \a *map on return contains a list of NoAtom entries such that going through the list, yields indices to the
+ * lines in \a *FileBuffer in a sorted manner of the Ion_Type?_? keywords. We assume that ConfigFileBuffer::CurrentLine
+ * points to first Ion_Type entry.
+ * \param *FileBuffer pointer to buffer structure
+ * \param NoAtoms of subsequent lines to look at
+ */
+void ConfigFileBuffer::MapIonTypesInBuffer(const int NoAtoms)
+{
+  map<const char *, int, IonTypeCompare> IonTypeLineMap;
+  if (!MappingAllocated) {
+    InitMapping();
+  }
+
+  // put all into hashed map
+  for (int i=0; i<NoAtoms; ++i) {
+    IonTypeLineMap.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i));
+  }
+
+  // fill map
+  int nr=0;
+  for (map<const char *, int, IonTypeCompare>::iterator runner = IonTypeLineMap.begin(); runner != IonTypeLineMap.end(); ++runner) {
+    if (CurrentLine+nr < NoLines)
+      LineMapping[CurrentLine+(nr++)] = runner->second;
+    else {
+      DoeLog(0) && (eLog()<< Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl);
+      performCriticalExit();
+    }
+  }
+}
Index: src/ConfigFileBuffer.hpp
===================================================================
--- src/ConfigFileBuffer.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/ConfigFileBuffer.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,53 @@
+/*
+ * ConfigFileBuffer.hpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#ifndef CONFIGFILEBUFFER_HPP_
+#define CONFIGFILEBUFFER_HPP_
+
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iostream>
+
+/****************************************** forward declarations *****************************/
+
+/********************************************** declarations *******************************/
+
+class ConfigFileBuffer {
+  public:
+    char **buffer;
+    int *LineMapping;
+    int CurrentLine;
+    int NoLines;
+
+    ConfigFileBuffer();
+    ConfigFileBuffer(const char * const filename);
+    ~ConfigFileBuffer();
+
+    void InitFileBuffer(const char * const filename);
+    void InitFileBuffer(std::istream *file);
+    void MapIonTypesInBuffer(const int NoAtoms);
+
+private:
+
+    /*
+     * Whether LineMapping is allocated or not.
+     */
+    bool MappingAllocated;
+
+    void InitMapping();
+    void RemoveMapping();
+    void RemoveBuffer();
+};
+
+
+#endif /* CONFIGFILEBUFFER_HPP_ */
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Makefile.am	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -88,9 +88,14 @@
   Parser/ChangeTracker.cpp \
   Parser/FormatParser.cpp \
+  Parser/MpqcParser.cpp \
+  Parser/PcpParser.cpp \
   Parser/TremoloParser.cpp \
   Parser/XyzParser.cpp
+
 PARSERHEADER = \
   Parser/ChangeTracker.hpp \
   Parser/FormatParser.hpp \
+  Parser/MpqcParser.hpp \
+  Parser/PcpParser.hpp \
   Parser/TremoloParser.hpp \
   Parser/XyzParser.hpp
@@ -153,4 +158,5 @@
   CommandLineParser.cpp \
   config.cpp \
+  ConfigFileBuffer.cpp \
   element.cpp \
   elements_db.cpp \
@@ -178,4 +184,5 @@
   tesselation.cpp \
   tesselationhelpers.cpp \
+  ThermoStatContainer.cpp \
   triangleintersectionlist.cpp \
   vector.cpp \
@@ -198,4 +205,5 @@
   CommandLineParser.hpp \
   config.hpp \
+  ConfigFileBuffer.hpp \
   defs.hpp \
   element.hpp \
@@ -220,4 +228,5 @@
   tesselation.hpp \
   tesselationhelpers.hpp \
+  ThermoStatContainer.hpp \
   triangleintersectionlist.hpp \
   verbose.hpp \
Index: src/Parser/ChangeTracker.cpp
===================================================================
--- src/Parser/ChangeTracker.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Parser/ChangeTracker.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -7,8 +7,7 @@
 
 #include "Helpers/MemDebug.hpp"
+#include "Parser/ChangeTracker.hpp"
+#include "Patterns/Singleton_impl.hpp"
 
-#include "ChangeTracker.hpp"
-
-ChangeTracker* ChangeTracker::instance = NULL;
 
 /**
@@ -27,26 +26,4 @@
 ChangeTracker::~ChangeTracker() {
   World::getInstance().signOff(this);
-}
-
-/**
- * Returns the change tracker instance.
- *
- * \return this
- */
-ChangeTracker* ChangeTracker::get() {
-  if (instance == NULL) {
-    instance = new ChangeTracker();
-  }
-
-  return instance;
-}
-
-/**
- * Destroys the change tracker instance. Be careful, the change tracker is a
- * singleton and destruction might lead to a loss of consistency.
- */
-void ChangeTracker::destroy() {
-  delete instance;
-  instance = NULL;
 }
 
@@ -76,2 +53,4 @@
   }
 }
+
+CONSTRUCT_SINGLETON(ChangeTracker)
Index: src/Parser/ChangeTracker.hpp
===================================================================
--- src/Parser/ChangeTracker.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Parser/ChangeTracker.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -18,5 +18,6 @@
  * changes to it.
  */
-class ChangeTracker : public Observable {
+class ChangeTracker : public Singleton<ChangeTracker>, public Observable {
+  friend class Singleton<ChangeTracker>;
 public:
   void saveStatus();
@@ -31,4 +32,6 @@
   bool isConsistent;
   static ChangeTracker* instance;
+
+  // private constructor and destructor due to singleton
   ChangeTracker();
   ~ChangeTracker();
Index: src/Parser/FormatParser.cpp
===================================================================
--- src/Parser/FormatParser.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Parser/FormatParser.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -19,5 +19,5 @@
   Observer("FormatParser")
 {
-  ChangeTracker::get()->signOn(this);
+  ChangeTracker::getInstance().signOn(this);
   saveStream = NULL;
 }
@@ -27,5 +27,5 @@
  */
 FormatParser::~FormatParser() {
-  ChangeTracker::get()->signOff(this);
+  ChangeTracker::getInstance().signOff(this);
 }
 
Index: src/Parser/MpqcParser.cpp
===================================================================
--- src/Parser/MpqcParser.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/Parser/MpqcParser.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,144 @@
+/*
+ * MpqcParser.cpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#include <iostream>
+
+#include "MpqcParser.hpp"
+
+#include "atom.hpp"
+#include "config.hpp"
+#include "element.hpp"
+#include "log.hpp"
+#include "periodentafel.hpp"
+#include "vector.hpp"
+#include "verbose.hpp"
+#include "World.hpp"
+
+
+/** Constructor of MpqcParser.
+ *
+ */
+MpqcParser::MpqcParser()
+{
+
+}
+
+/** Destructor of MpqcParser.
+ *
+ */
+MpqcParser::~MpqcParser() {
+
+}
+
+/** Load an MPQC config file into the World.
+ * \param *file input stream
+ */
+void MpqcParser::load(istream *file)
+{
+  DoeLog(0) && (Log() << Verbose(0) << "Not yet implemented" << endl) ;
+}
+
+void MpqcParser::save(ostream *file)
+{
+  if (HessianPresent)
+    saveHessian(file);
+  else
+    saveSimple(file);
+}
+
+/** Saves all atoms and data into a MPQC config file without hessian.
+ * \param *file output stream
+ */
+void MpqcParser::saveSimple(ostream *file)
+{
+  int AtomNo = 0;
+  Vector center;
+  vector<atom *> allatoms = World::getInstance().getAllAtoms();
+
+  // first without hessian
+  if (file->fail()) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc output file." << endl);
+  } else {
+    *file << "% Created by MoleCuilder" << endl;
+    *file << "mpqc: (" << endl;
+    *file << "\tsavestate = no" << endl;
+    *file << "\tdo_gradient = yes" << endl;
+    *file << "\tmole<MBPT2>: (" << endl;
+    *file << "\t\tmaxiter = 200" << endl;
+    *file << "\t\tbasis = $:basis" << endl;
+    *file << "\t\tmolecule = $:molecule" << endl;
+    *file << "\t\treference<CLHF>: (" << endl;
+    *file << "\t\t\tbasis = $:basis" << endl;
+    *file << "\t\t\tmolecule = $:molecule" << endl;
+    *file << "\t\t)" << endl;
+    *file << "\t)" << endl;
+    *file << ")" << endl;
+    *file << "molecule<Molecule>: (" << endl;
+    *file << "\tunit = " << (World::getInstance().getConfig()->GetIsAngstroem() ? "angstrom" : "bohr" ) << endl;
+    *file << "\t{ atoms geometry } = {" << endl;
+    // output of atoms
+    for (vector<atom *>::iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) {
+      (*AtomRunner)->OutputMPQCLine(file, &center, &AtomNo);
+    }
+    *file << "\t}" << endl;
+    *file << ")" << endl;
+    *file << "basis<GaussianBasisSet>: (" << endl;
+    *file << "\tname = \"" << World::getInstance().getConfig()->basis << "\"" << endl;
+    *file << "\tmolecule = $:molecule" << endl;
+    *file << ")" << endl;
+  }
+}
+
+/** Saves all atoms and data into a MPQC config file with hessian.
+ * \param *file output stream
+ */
+void MpqcParser::saveHessian(ostream *file)
+{
+  int AtomNo = 0;
+  Vector center;
+  vector<atom *> allatoms = World::getInstance().getAllAtoms();
+
+  // with hessian
+  if (file->fail()) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc output file." << endl);
+  } else {
+    *file << "% Created by MoleCuilder" << endl;
+    *file << "mpqc: (" << endl;
+    *file << "\tsavestate = no" << endl;
+    *file << "\tdo_gradient = yes" << endl;
+    *file << "\tmole<CLHF>: (" << endl;
+    *file << "\t\tmaxiter = 200" << endl;
+    *file << "\t\tbasis = $:basis" << endl;
+    *file << "\t\tmolecule = $:molecule" << endl;
+    *file << "\t)" << endl;
+    *file << "\tfreq<MolecularFrequencies>: (" << endl;
+    *file << "\t\tmolecule=$:molecule" << endl;
+    *file << "\t)" << endl;
+    *file << ")" << endl;
+    *file << "molecule<Molecule>: (" << endl;
+    *file << "\tunit = " << (World::getInstance().getConfig()->GetIsAngstroem() ? "angstrom" : "bohr" ) << endl;
+    *file << "\t{ atoms geometry } = {" << endl;
+    // output of atoms
+    for (vector<atom *>::iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) {
+      (*AtomRunner)->OutputMPQCLine(file, &center, &AtomNo);
+    }
+    *file << "\t}" << endl;
+    *file << ")" << endl;
+    *file << "basis<GaussianBasisSet>: (" << endl;
+    *file << "\tname = \"" << World::getInstance().getConfig()->basis << "\"" << endl;
+    *file << "\tmolecule = $:molecule" << endl;
+    *file << ")" << endl;
+  }
+}
+
+/** Sets whether hessian is desired or not
+ * \param hessian statement
+ */
+void MpqcParser::setHessian(bool hessian)
+{
+  HessianPresent = hessian;
+}
Index: src/Parser/MpqcParser.hpp
===================================================================
--- src/Parser/MpqcParser.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/Parser/MpqcParser.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,35 @@
+/*
+ * MpqcParser.hpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#ifndef MPQCPARSER_HPP_
+#define MPQCPARSER_HPP_
+
+#include <iostream>
+#include "FormatParser.hpp"
+
+/**
+ * Loads a MPQC config file into the World and saves the World as a tremolo file.
+ */
+class MpqcParser : public FormatParser
+{
+public:
+  MpqcParser();
+  ~MpqcParser();
+  void load(std::istream* file);
+  void save(std::ostream* file);
+
+  void setHessian(bool hessian);
+
+private:
+  // whether output with computing hessian is desired.
+  bool HessianPresent;
+
+  void saveSimple(std::ostream *file);
+  void saveHessian(std::ostream *file);
+};
+
+#endif /* MPQCPARSER_HPP_ */
Index: src/Parser/PcpParser.cpp
===================================================================
--- src/Parser/PcpParser.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/Parser/PcpParser.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,648 @@
+/*
+ * PcpParser.cpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#include <iostream>
+
+#include "atom.hpp"
+#include "config.hpp"
+#include "ConfigFileBuffer.hpp"
+#include "element.hpp"
+#include "Helpers/Assert.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()
+{
+  Parallelization.ProcPEGamma = 1;
+  Parallelization.ProcPEPsi = 1;
+
+  Paths.databasepath = NULL;
+  Paths.configpath = NULL;
+  Paths.configname = NULL;
+  Paths.mainname = NULL;
+  Paths.defaultpath = NULL;
+  Paths.pseudopotpath = NULL;
+
+  Switches.DoConstrainedMD = 0;
+  Switches.DoOutVis = 0;
+  Switches.DoOutMes = 1;
+  Switches.DoOutNICS = 0;
+  Switches.DoOutOrbitals = 0;
+  Switches.DoOutCurrent = 0;
+  Switches.DoFullCurrent = 0;
+  Switches.DoPerturbation = 0;
+  Switches.DoWannier = 0;
+
+  LocalizedOrbitals.CommonWannier = 0;
+  LocalizedOrbitals.SawtoothStart = 0.01;
+  LocalizedOrbitals.VectorPlane = 0;
+  LocalizedOrbitals.VectorCut = 0;
+  LocalizedOrbitals.UseAddGramSch = 1;
+  LocalizedOrbitals.Seed = 1;
+  LocalizedOrbitals.EpsWannier = 1e-7;
+
+  StepCounts.MaxMinStopStep = 0;
+  StepCounts.InitMaxMinStopStep = 0;
+  StepCounts.OutVisStep = 10;
+  StepCounts.OutSrcStep = 5;
+  StepCounts.MaxPsiStep = 3;
+  StepCounts.MaxOuterStep = 0;
+  StepCounts.MaxMinStep = 100;
+  StepCounts.RelEpsTotalEnergy = 1e-07;
+  StepCounts.RelEpsKineticEnergy = 1e-05;
+  StepCounts.MaxMinGapStopStep = 0;
+  StepCounts.MaxInitMinStep = 100;
+  StepCounts.InitRelEpsTotalEnergy = 1e-05;
+  StepCounts.InitRelEpsKineticEnergy = 0.0001;
+  StepCounts.InitMaxMinGapStopStep = 1;
+
+  PlaneWaveSpecifics.PsiType = 0;
+  PlaneWaveSpecifics.MaxPsiDouble = 0;
+  PlaneWaveSpecifics.PsiMaxNoUp = 0;
+  PlaneWaveSpecifics.PsiMaxNoDown = 0;
+  PlaneWaveSpecifics.ECut = 128;
+  PlaneWaveSpecifics.MaxLevel = 5;
+  PlaneWaveSpecifics.RiemannTensor = 0;
+  PlaneWaveSpecifics.LevRFactor = 0;
+  PlaneWaveSpecifics.RiemannLevel = 0;
+  PlaneWaveSpecifics.Lev0Factor = 2;
+  PlaneWaveSpecifics.RTActualUse = 0;
+  PlaneWaveSpecifics.AddPsis = 0;
+  PlaneWaveSpecifics.RCut = 20;
+  PlaneWaveSpecifics.PsiType = 0;
+
+  FastParsing = false;
+
+  Deltat = 0.01;
+  IsAngstroem = 1;
+  RelativeCoord = 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 * 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;
+  int nr = 0;
+  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;
+    *file << "Ion_Type" << ZtoIndexMap[(*AtomRunner)->type->Z] << "_" << ZtoCountMap[Inserter.first->second] << "\t"  << fixed << setprecision(9) << showpoint;
+    *file << (*AtomRunner)->x[0] << "\t" << (*AtomRunner)->x[1] << "\t" << (*AtomRunner)->x[2];
+    *file << "\t" << (*AtomRunner)->FixedIon;
+    if ((*AtomRunner)->v.Norm() > MYEPSILON)
+      *file << "\t" << scientific << setprecision(6) << (*AtomRunner)->v[0] << "\t" << (*AtomRunner)->v[1] << "\t" << (*AtomRunner)->v[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)) {
+    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);
+};
+
+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.configpath != NULL) && (b.Paths.configpath != NULL))
+    ASSERT(strcmp(Paths.configpath, b.Paths.configpath), "PcpParser ==: configpath 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;
+}
Index: src/Parser/PcpParser.hpp
===================================================================
--- src/Parser/PcpParser.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/Parser/PcpParser.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,129 @@
+/*
+ * PcpParser.hpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#ifndef PCPPARSER_HPP_
+#define PCPPARSER_HPP_
+
+#include <iostream>
+#include "FormatParser.hpp"
+
+/**
+ * Parser for PCP config files.
+ */
+class PcpParser : public FormatParser {
+public:
+  PcpParser();
+  ~PcpParser();
+  void load(std::istream* file);
+  void save(std::ostream* file);
+
+  bool operator==(const PcpParser& b) const;
+
+private:
+
+  void ParseThermostats(class ConfigFileBuffer * const fb);
+  void OutputAtoms(ostream *file, vector<atom *> &allatoms, map<int, int> &ZtoIndexMap);
+  void OutputElements(ostream *file, vector<atom *> &allatoms, map<int, int> &ZtoIndexMap);
+
+  struct StructParallelization {
+    int ProcPEGamma;
+    int ProcPEPsi;
+  } Parallelization;
+
+  /*
+   * Contains all the paths and names
+   */
+  struct StructPaths {
+    char *databasepath;
+    char *configpath;
+    char *configname;
+    char *mainname;
+    char *defaultpath;
+    char *pseudopotpath;
+  } Paths;
+
+  /*
+   * Contains all Do/Don't switches
+   */
+  struct StructSwitches {
+    int DoConstrainedMD;
+    int DoOutVis;
+    int DoOutMes;
+    int DoOutNICS;
+    int DoOutOrbitals;
+    int DoOutCurrent;
+    int DoFullCurrent;
+    int DoPerturbation;
+    int DoWannier;
+  } Switches;
+
+  /*
+   * Contains parameters regarding localization of orbitals or magnetic perturbation
+   */
+  struct StructLocalizedOrbitals {
+    int CommonWannier;
+    double SawtoothStart;
+    int VectorPlane;
+    double VectorCut;
+    int UseAddGramSch;
+    int Seed;
+    double EpsWannier;
+  } LocalizedOrbitals;
+
+  /*
+   * Contains all step count and other epsilon threshold parameters
+   */
+  struct StructStepCounts {
+    int MaxMinStopStep;
+    int InitMaxMinStopStep;
+
+    int OutVisStep;
+    int OutSrcStep;
+    int MaxPsiStep;
+
+    int MaxOuterStep;
+
+    int MaxMinStep;
+    double RelEpsTotalEnergy;
+    double RelEpsKineticEnergy;
+    int MaxMinGapStopStep;
+    int MaxInitMinStep;
+    double InitRelEpsTotalEnergy;
+    double InitRelEpsKineticEnergy;
+    int InitMaxMinGapStopStep;
+  } StepCounts;
+
+  /*
+   * Contains all parameters specific to the plane wave basis set
+   */
+  struct StructPlaneWaveSpecifics {
+    int PsiType;
+    int MaxPsiDouble;
+    int PsiMaxNoUp;
+    int PsiMaxNoDown;
+    double ECut;
+    int MaxLevel;
+    int RiemannTensor;
+    int LevRFactor;
+    int RiemannLevel;
+    int Lev0Factor;
+    int RTActualUse;
+    int AddPsis;
+    double RCut;
+  } PlaneWaveSpecifics;
+
+  bool FastParsing;
+
+  double Deltat;
+  int IsAngstroem;
+  int RelativeCoord;
+  int StructOpt;
+  int MaxTypes;
+  string basis;
+};
+
+#endif /* PCPPARSER_HPP_ */
Index: src/Parser/TremoloParser.hpp
===================================================================
--- src/Parser/TremoloParser.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Parser/TremoloParser.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -10,5 +10,5 @@
 
 #include <string>
-#include "FormatParser.hpp"
+#include "Parser/FormatParser.hpp"
 
 /**
Index: src/Parser/XyzParser.cpp
===================================================================
--- src/Parser/XyzParser.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Parser/XyzParser.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -27,5 +27,4 @@
  */
 XyzParser::~XyzParser() {
-  delete(&comment);
 }
 
Index: src/Parser/XyzParser.hpp
===================================================================
--- src/Parser/XyzParser.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/Parser/XyzParser.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -10,5 +10,5 @@
 
 #include <string>
-#include "FormatParser.hpp"
+#include "Parser/FormatParser.hpp"
 
 /**
Index: src/ThermoStatContainer.cpp
===================================================================
--- src/ThermoStatContainer.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/ThermoStatContainer.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,52 @@
+/*
+ * ThermoStatContainer.cpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#include <cstring>
+
+#include "ConfigFileBuffer.hpp"
+#include "config.hpp"
+#include "log.hpp"
+#include "ThermoStatContainer.hpp"
+#include "verbose.hpp"
+
+/** Constructor for class ThermoStatContainer.
+ *
+ */
+ThermoStatContainer::ThermoStatContainer() : Thermostat(4), ThermostatImplemented(NULL),
+ThermostatNames(NULL), TempFrequency(2.5), alpha(0.), HooverMass(0.), TargetTemp(0.00095004455), ScaleTempStep(25)
+{
+  ThermostatImplemented = new int[MaxThermostats];
+  ThermostatNames = new char *[MaxThermostats];
+  for (int j=0;j<MaxThermostats;j++)
+    ThermostatNames[j] = new char[12];
+
+  strcpy(ThermostatNames[0],"None");
+  ThermostatImplemented[0] = 1;
+  strcpy(ThermostatNames[1],"Woodcock");
+  ThermostatImplemented[1] = 1;
+  strcpy(ThermostatNames[2],"Gaussian");
+  ThermostatImplemented[2] = 1;
+  strcpy(ThermostatNames[3],"Langevin");
+  ThermostatImplemented[3] = 1;
+  strcpy(ThermostatNames[4],"Berendsen");
+  ThermostatImplemented[4] = 1;
+  strcpy(ThermostatNames[5],"NoseHoover");
+  ThermostatImplemented[5] = 1;
+}
+
+/** Destructor for Class ThermoStatContainer.
+ *
+ */
+ThermoStatContainer::~ThermoStatContainer()
+{
+  delete[](ThermostatImplemented);
+  for (int j=0;j<MaxThermostats;j++)
+    delete[](ThermostatNames[j]);
+  delete[](ThermostatNames);
+}
+
+
Index: src/ThermoStatContainer.hpp
===================================================================
--- src/ThermoStatContainer.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
+++ src/ThermoStatContainer.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -0,0 +1,38 @@
+/*
+ * ThermoStatContainer.hpp
+ *
+ *  Created on: 12.06.2010
+ *      Author: heber
+ */
+
+#ifndef THERMOSTATCONTAINER_HPP_
+#define THERMOSTATCONTAINER_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+class ConfigFileBuffer;
+
+#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
+enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
+
+class ThermoStatContainer {
+public:
+  ThermoStatContainer();
+  ~ThermoStatContainer();
+
+  void ParseThermostats(class ConfigFileBuffer * const fb);
+
+  int Thermostat;
+  int *ThermostatImplemented;
+  char **ThermostatNames;
+  double TempFrequency;
+  double alpha;
+  double HooverMass;
+  double TargetTemp;
+  int ScaleTempStep;
+};
+
+#endif /* THERMOSTATCONTAINER_HPP_ */
Index: src/UIElements/CommandLineUI/CommandLineWindow.cpp
===================================================================
--- src/UIElements/CommandLineUI/CommandLineWindow.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/UIElements/CommandLineUI/CommandLineWindow.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -83,7 +83,11 @@
   // go through all possible actions
   for (std::list<std::string>::iterator CommandRunner = CommandLineParser::getInstance().SequenceOfActions.begin(); CommandRunner != CommandLineParser::getInstance().SequenceOfActions.end(); ++CommandRunner) {
-    cout << "Checking presence of " << *CommandRunner << endl;
-    if (ActionRegistry::getInstance().isActionByNamePresent(*CommandRunner))
+    cout << "Checking presence of " << *CommandRunner << ": ";
+    if (ActionRegistry::getInstance().isActionByNamePresent(*CommandRunner)) {
+      cout << "calling " << *CommandRunner << endl;
       ActionRegistry::getInstance().getActionByName(*CommandRunner)->call();
+    } else {
+      cout << "absent." << endl;
+    }
   }
 }
Index: src/World.cpp
===================================================================
--- src/World.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/World.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -14,4 +14,5 @@
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "ThermoStatContainer.hpp"
 #include "Descriptors/AtomDescriptor.hpp"
 #include "Descriptors/AtomDescriptor_impl.hpp"
@@ -89,4 +90,10 @@
   defaultName = name;
 };
+
+class ThermoStatContainer * World::getThermostats()
+{
+  return Thermostats;
+}
+
 
 int World::getExitFlag() {
@@ -281,4 +288,5 @@
     periode(new periodentafel),
     configuration(new config),
+    Thermostats(new ThermoStatContainer),
     ExitFlag(0),
     atoms(),
@@ -306,4 +314,5 @@
   delete periode;
   delete configuration;
+  delete Thermostats;
   MoleculeSet::iterator molIter;
   for(molIter=molecules.begin();molIter!=molecules.end();++molIter){
Index: src/World.hpp
===================================================================
--- src/World.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/World.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -30,16 +30,16 @@
 
 // forward declarations
-class config;
-class periodentafel;
-class MoleculeListClass;
 class atom;
-class molecule;
 class AtomDescriptor;
 class AtomDescriptor_impl;
+template<typename T> class AtomsCalculation;
+class config;
+class ManipulateAtomsProcess;
+class molecule;
 class MoleculeDescriptor;
 class MoleculeDescriptor_impl;
-class ManipulateAtomsProcess;
-template<typename T>
-class AtomsCalculation;
+class MoleculeListClass;
+class periodentafel;
+class ThermoStatContainer;
 
 /****************************************** forward declarations *****************************/
@@ -142,4 +142,9 @@
   void setDefaultName(std::string name);
 
+  /**
+   * get pointer to World's ThermoStatContainer
+   */
+  ThermoStatContainer * getThermostats();
+
   /*
    * get the ExitFlag
@@ -254,4 +259,5 @@
   static double *cell_size;
   std::string defaultName;
+  class ThermoStatContainer *Thermostats;
   int ExitFlag;
 public:
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/atom.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -159,5 +159,5 @@
   * \return true - \a *out present, false - \a *out is NULL
  */
-bool atom::OutputArrayIndexed(ofstream * const out, const int *ElementNo, int *AtomNo, const char *comment) const
+bool atom::OutputArrayIndexed(ostream * const out, const int *ElementNo, int *AtomNo, const char *comment) const
 {
   AtomNo[type->Z]++;  // increment number
@@ -236,5 +236,5 @@
  * \param *AtomNo pointer to atom counter that is increased by one
  */
-void atom::OutputMPQCLine(ofstream * const out, const Vector *center, int *AtomNo = NULL) const
+void atom::OutputMPQCLine(ostream * const out, const Vector *center, int *AtomNo = NULL) const
 {
   *out << "\t\t" << type->symbol << " [ " << x[0]-center->at(0) << "\t" << x[1]-center->at(1) << "\t" << x[2]-center->at(2) << " ]" << endl;
Index: src/atom.hpp
===================================================================
--- src/atom.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/atom.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -51,9 +51,9 @@
 
   bool OutputIndexed(ofstream * const out, const int ElementNo, const int AtomNo, const char *comment = NULL) const;
-  bool OutputArrayIndexed(ofstream * const out, const int *ElementNo, int *AtomNo, const char *comment = NULL) const;
+  bool OutputArrayIndexed(ostream * const out, const int *ElementNo, int *AtomNo, const char *comment = NULL) const;
   bool OutputXYZLine(ofstream *out) const;
   bool OutputTrajectory(ofstream * const out, const int *ElementNo, int *AtomNo, const int step) const;
   bool OutputTrajectoryXYZ(ofstream * const out, const int step) const;
-  void OutputMPQCLine(ofstream * const out, const Vector *center, int *AtomNo) const;
+  void OutputMPQCLine(ostream * const out, const Vector *center, int *AtomNo) const;
 
   void InitComponentNr();
Index: src/atom_trajectoryparticle.cpp
===================================================================
--- src/atom_trajectoryparticle.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/atom_trajectoryparticle.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -15,4 +15,5 @@
 #include "log.hpp"
 #include "parser.hpp"
+#include "ThermoStatContainer.hpp"
 #include "verbose.hpp"
 
@@ -197,9 +198,9 @@
 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
 {
-  double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+  double sigma  = sqrt(configuration->Thermostats->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
   Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     // throw a dice to determine whether it gets hit by a heat bath particle
-    if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
+    if (((((rand()/(double)RAND_MAX))*configuration->Thermostats->TempFrequency) < 1.)) {
       DoLog(3) && (Log() << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ");
       // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
@@ -225,5 +226,5 @@
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
-      U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
+      U[d] *= sqrt(1+(configuration->Deltat/configuration->Thermostats->TempFrequency)*(ScaleTempFactor-1));
       *ekin += 0.5*type->mass * U[d]*U[d];
     }
@@ -255,5 +256,5 @@
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
-        U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
+        U[d] += configuration->Deltat/type->mass * (configuration->Thermostats->alpha * (U[d] * type->mass));
         *ekin += (0.5*type->mass) * U[d]*U[d];
       }
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/builder.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -79,4 +79,5 @@
 #include "UIElements/Dialog.hpp"
 #include "Menu/ActionMenuItem.hpp"
+#include "Parser/ChangeTracker.hpp"
 #include "Actions/ActionRegistry.hpp"
 #include "Actions/ActionHistory.hpp"
@@ -1691,11 +1692,8 @@
               } else {
                 SaveFlag = true;
-                DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl);
-                if (!mol->AddXYZFile(argv[argptr]))
-                  DoLog(2) && (Log() << Verbose(2) << "File not found." << endl);
-                else {
-                  DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl);
-                  configPresent = present;
-                }
+                configPresent = present;
+                ArgcList.insert(argptr-1);
+                ArgcList.insert(argptr);
+                argptr+=1;
               }
               break;
@@ -2181,4 +2179,5 @@
 
 void cleanUp(){
+  ChangeTracker::purgeInstance();
   World::purgeInstance();
   logger::purgeInstance();
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/config.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -10,8 +10,9 @@
 #include <cstring>
 
-#include "World.hpp"
 #include "atom.hpp"
 #include "bond.hpp"
+#include "bondgraph.hpp"
 #include "config.hpp"
+#include "ConfigFileBuffer.hpp"
 #include "element.hpp"
 #include "helpers.hpp"
@@ -23,150 +24,6 @@
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "ThermoStatContainer.hpp"
 #include "World.hpp"
-
-/******************************** Functions for class ConfigFileBuffer **********************/
-
-/** Structure containing compare function for Ion_Type sorting.
- */
-struct IonTypeCompare {
-  bool operator()(const char* s1, const char *s2) const {
-    char number1[8];
-    char number2[8];
-    const char *dummy1, *dummy2;
-    //Log() << Verbose(0) << s1 << "  " << s2 << endl;
-    dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
-    dummy2 = strchr(dummy1, '_');
-    strncpy(number1, dummy1, dummy2-dummy1); // copy the number
-    number1[dummy2-dummy1]='\0';
-    dummy1 = strchr(s2, '_')+sizeof(char)*5;  // go just after "Ion_Type"
-    dummy2 = strchr(dummy1, '_');
-    strncpy(number2, dummy1, dummy2-dummy1); // copy the number
-    number2[dummy2-dummy1]='\0';
-    if (atoi(number1) != atoi(number2))
-      return (atoi(number1) < atoi(number2));
-    else {
-      dummy1 = strchr(s1, '_')+sizeof(char);
-      dummy1 = strchr(dummy1, '_')+sizeof(char);
-      dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
-      strncpy(number1, dummy1, dummy2-dummy1); // copy the number
-      number1[dummy2-dummy1]='\0';
-      dummy1 = strchr(s2, '_')+sizeof(char);
-      dummy1 = strchr(dummy1, '_')+sizeof(char);
-      dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
-      strncpy(number2, dummy1, dummy2-dummy1); // copy the number
-      number2[dummy2-dummy1]='\0';
-      return (atoi(number1) < atoi(number2));
-    }
-  }
-};
-
-/** Constructor for ConfigFileBuffer class.
- */
-ConfigFileBuffer::ConfigFileBuffer() : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
-{
-};
-
-/** Constructor for ConfigFileBuffer class with filename to be parsed.
- * \param *filename file name
- */
-ConfigFileBuffer::ConfigFileBuffer(const char * const filename) : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
-{
-  ifstream *file = NULL;
-  char line[MAXSTRINGSIZE];
-
-  // prescan number of lines
-  file= new ifstream(filename);
-  if (file == NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
-    return;
-  }
-  NoLines = 0; // we're overcounting by one
-  long file_position = file->tellg(); // mark current position
-  do {
-    file->getline(line, 256);
-    NoLines++;
-  } while (!file->eof());
-  file->clear();
-  file->seekg(file_position, ios::beg);
-  DoLog(1) && (Log() << Verbose(1) << NoLines-1 << " lines were recognized." << endl);
-
-  // allocate buffer's 1st dimension
-  if (buffer != NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "FileBuffer->buffer is not NULL!" << endl);
-    return;
-  } else
-    buffer = new char *[NoLines];
-
-  // scan each line and put into buffer
-  int lines=0;
-  int i;
-  do {
-    buffer[lines] = new char[MAXSTRINGSIZE];
-    file->getline(buffer[lines], MAXSTRINGSIZE-1);
-    i = strlen(buffer[lines]);
-    buffer[lines][i] = '\n';
-    buffer[lines][i+1] = '\0';
-    lines++;
-  } while((!file->eof()) && (lines < NoLines));
-  DoLog(1) && (Log() << Verbose(1) << lines-1 << " lines were read into the buffer." << endl);
-
-  // close and exit
-  file->close();
-  file->clear();
-  delete(file);
-}
-
-/** Destructor for ConfigFileBuffer class.
- */
-ConfigFileBuffer::~ConfigFileBuffer()
-{
-  for(int i=0;i<NoLines;++i)
-    delete[](buffer[i]);
-  delete[](buffer);
-  delete[](LineMapping);
-}
-
-
-/** Create trivial mapping.
- */
-void ConfigFileBuffer::InitMapping()
-{
-  LineMapping = new int[NoLines];
-  for (int i=0;i<NoLines;i++)
-    LineMapping[i] = i;
-}
-
-/** Creates a mapping for the \a *FileBuffer's lines containing the Ion_Type keyword such that they are sorted.
- * \a *map on return contains a list of NoAtom entries such that going through the list, yields indices to the
- * lines in \a *FileBuffer in a sorted manner of the Ion_Type?_? keywords. We assume that ConfigFileBuffer::CurrentLine
- * points to first Ion_Type entry.
- * \param *FileBuffer pointer to buffer structure
- * \param NoAtoms of subsequent lines to look at
- */
-void ConfigFileBuffer::MapIonTypesInBuffer(const int NoAtoms)
-{
-  map<const char *, int, IonTypeCompare> IonTypeLineMap;
-  if (LineMapping == NULL) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "map pointer is NULL: " << LineMapping << endl);
-    performCriticalExit();
-    return;
-  }
-
-  // put all into hashed map
-  for (int i=0; i<NoAtoms; ++i) {
-    IonTypeLineMap.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i));
-  }
-
-  // fill map
-  int nr=0;
-  for (map<const char *, int, IonTypeCompare>::iterator runner = IonTypeLineMap.begin(); runner != IonTypeLineMap.end(); ++runner) {
-    if (CurrentLine+nr < NoLines)
-      LineMapping[CurrentLine+(nr++)] = runner->second;
-    else {
-      DoeLog(0) && (eLog()<< Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl);
-      performCriticalExit();
-    }
-  }
-}
 
 /************************************* Functions for class config ***************************/
@@ -174,7 +31,6 @@
 /** Constructor for config file class.
  */
-config::config() : BG(NULL), PsiType(0), MaxPsiDouble(0), PsiMaxNoUp(0), PsiMaxNoDown(0), MaxMinStopStep(1), InitMaxMinStopStep(1), ProcPEGamma(8), ProcPEPsi(1), configpath(NULL),
-    configname(NULL), FastParsing(false), Deltat(0.01), basis(""), databasepath(NULL), DoConstrainedMD(0), MaxOuterStep(0), Thermostat(4), ThermostatImplemented(NULL),
-    ThermostatNames(NULL), TempFrequency(2.5), alpha(0.), HooverMass(0.), TargetTemp(0.00095004455), ScaleTempStep(25),  mainname(NULL), defaultpath(NULL), pseudopotpath(NULL),
+config::config() : BG(NULL), Thermostats(0), PsiType(0), MaxPsiDouble(0), PsiMaxNoUp(0), PsiMaxNoDown(0), MaxMinStopStep(1), InitMaxMinStopStep(1), ProcPEGamma(8), ProcPEPsi(1), configpath(NULL),
+    configname(NULL), FastParsing(false), Deltat(0.01), basis(""), databasepath(NULL), DoConstrainedMD(0), MaxOuterStep(0), mainname(NULL), defaultpath(NULL), pseudopotpath(NULL),
     DoOutVis(0), DoOutMes(1), DoOutNICS(0), DoOutOrbitals(0), DoOutCurrent(0), DoFullCurrent(0), DoPerturbation(0), DoWannier(0), CommonWannier(0), SawtoothStart(0.01),
     VectorPlane(0), VectorCut(0.), UseAddGramSch(1), Seed(1), OutVisStep(10), OutSrcStep(5), MaxPsiStep(0), EpsWannier(1e-7), MaxMinStep(100), RelEpsTotalEnergy(1e-7),
@@ -188,4 +44,5 @@
   configpath = new char[MAXSTRINGSIZE];
   configname = new char[MAXSTRINGSIZE];
+  Thermostats = new ThermoStatContainer();
   strcpy(mainname,"pcp");
   strcpy(defaultpath,"not specified");
@@ -194,6 +51,4 @@
   configname[0]='\0';
   basis = "3-21G";
-
-  InitThermostats();
 };
 
@@ -208,111 +63,10 @@
   delete[](configpath);
   delete[](configname);
-  delete[](ThermostatImplemented);
-  for (int j=0;j<MaxThermostats;j++)
-    delete[](ThermostatNames[j]);
-  delete[](ThermostatNames);
+  if (Thermostats != NULL)
+    delete(Thermostats);
 
   if (BG != NULL)
     delete(BG);
 };
-
-/** Initialises variables in class config for Thermostats.
- */
-void config::InitThermostats()
-{
-  ThermostatImplemented = new int[MaxThermostats];
-  ThermostatNames = new char *[MaxThermostats];
-  for (int j=0;j<MaxThermostats;j++)
-    ThermostatNames[j] = new char[12];
-
-  strcpy(ThermostatNames[0],"None");
-  ThermostatImplemented[0] = 1;
-  strcpy(ThermostatNames[1],"Woodcock");
-  ThermostatImplemented[1] = 1;
-  strcpy(ThermostatNames[2],"Gaussian");
-  ThermostatImplemented[2] = 1;
-  strcpy(ThermostatNames[3],"Langevin");
-  ThermostatImplemented[3] = 1;
-  strcpy(ThermostatNames[4],"Berendsen");
-  ThermostatImplemented[4] = 1;
-  strcpy(ThermostatNames[5],"NoseHoover");
-  ThermostatImplemented[5] = 1;
-};
-
-/** Readin of Thermostat related values from parameter file.
- * \param *fb file buffer containing the config file
- */
-void config::ParseThermostats(class ConfigFileBuffer * const fb)
-{
-  char * const thermo = new char[12];
-  const int verbose = 0;
-
-  // read desired Thermostat from file along with needed additional parameters
-  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
-    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
-      if (ThermostatImplemented[0] == 1) {
-        Thermostat = None;
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
-        Thermostat = None;
-      }
-    } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
-      if (ThermostatImplemented[1] == 1) {
-        Thermostat = Woodcock;
-        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
-        Thermostat = None;
-      }
-    } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
-      if (ThermostatImplemented[2] == 1) {
-        Thermostat = Gaussian;
-        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
-        Thermostat = None;
-      }
-    } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
-      if (ThermostatImplemented[3] == 1) {
-        Thermostat = Langevin;
-        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
-        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
-          DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl);
-        } else {
-          alpha = 1.;
-        }
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
-        Thermostat = None;
-      }
-    } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
-      if (ThermostatImplemented[4] == 1) {
-        Thermostat = Berendsen;
-        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
-        Thermostat = None;
-      }
-    } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
-      if (ThermostatImplemented[5] == 1) {
-        Thermostat = NoseHoover;
-        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
-        alpha = 0.;
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
-        Thermostat = None;
-      }
-    } else {
-      DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
-      Thermostat = None;
-    }
-  } else {
-    if ((MaxOuterStep > 0) && (TargetTemp != 0)) 
-      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
-    Thermostat = None;
-  }
-  delete[](thermo);
-};
-
 
 /** Displays menu for editing each entry of the config file.
@@ -656,19 +410,4 @@
 };
 
-/** Initializes ConfigFileBuffer from a file.
- * \param *file input file stream being the opened config file
- * \param *FileBuffer pointer to FileBuffer on return, should point to NULL
- */
-void PrepareFileBuffer(const char * const filename, struct ConfigFileBuffer *&FileBuffer)
-{
-  if (FileBuffer != NULL) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "deleting present FileBuffer in PrepareFileBuffer()." << endl);
-    delete(FileBuffer);
-  }
-  FileBuffer = new ConfigFileBuffer(filename);
-
-  FileBuffer->InitMapping();
-};
-
 /** Loads a molecule from a ConfigFileBuffer.
  * \param *mol molecule to load
@@ -867,6 +606,5 @@
 
   // ParseParameterFile
-  struct ConfigFileBuffer *FileBuffer = NULL;
-  PrepareFileBuffer(filename,FileBuffer);
+  class ConfigFileBuffer *FileBuffer = new ConfigFileBuffer(filename);
 
   /* Oeffne Hauptparameterdatei */
@@ -877,4 +615,5 @@
   int verbose = 0;
   
+  //TODO: This is actually sensible?: if (MaxOuterStep > 0)
   ParseThermostats(FileBuffer);
   
@@ -941,5 +680,5 @@
   ParseForParameter(verbose,FileBuffer,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
   ParseForParameter(verbose,FileBuffer,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
-  ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
+  ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(Thermostats->TargetTemp), 1, optional);
   //ParseForParameter(verbose,FileBuffer,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
   if (!ParseForParameter(verbose,FileBuffer,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
@@ -1150,6 +889,6 @@
   ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
   ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
-  ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
-  ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
+  ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(Thermostats->TargetTemp), 1, optional);
+  ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(Thermostats->ScaleTempStep), 1, optional);
   config::EpsWannier = 1e-8;
 
@@ -1339,23 +1078,23 @@
     *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
     *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
-    *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
-    switch(Thermostat) {
+    *output << "Thermostat\t" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
+    switch(Thermostats->Thermostat) {
       default:
       case None:
         break;
       case Woodcock:
-        *output << ScaleTempStep;
+        *output << Thermostats->ScaleTempStep;
         break;
       case Gaussian:
-        *output << ScaleTempStep;
+        *output << Thermostats->ScaleTempStep;
         break;
       case Langevin:
-        *output << TempFrequency << "\t" << alpha;
+        *output << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
         break;
       case Berendsen:
-        *output << TempFrequency;
+        *output << Thermostats->TempFrequency;
         break;
       case NoseHoover:
-        *output << HooverMass;
+        *output << Thermostats->HooverMass;
         break;
     };
@@ -1372,5 +1111,5 @@
     *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
     *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
-    *output << "TargetTemp\t" << config::TargetTemp << "\t# Target temperature" << endl;
+    *output << "TargetTemp\t" << Thermostats->TargetTemp << "\t# Target temperature" << endl;
     *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
     *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
@@ -1483,5 +1222,5 @@
     // output of atoms
     AtomNo = 0;
-    mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
+    mol->ActOnAllAtoms( &atom::OutputMPQCLine, (ostream * const) output, (const Vector *)center, &AtomNo );
     delete(center);
     *output << "\t}" << endl;
@@ -1525,5 +1264,5 @@
     // output of atoms
     AtomNo = 0;
-    mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
+    mol->ActOnAllAtoms( &atom::OutputMPQCLine, (ostream * const) output, (const Vector *)center, &AtomNo );
     delete(center);
     *output << "\t}" << endl;
@@ -2346,2 +2085,78 @@
   return (found); // true if found, false if not
 }
+
+/** Reading of Thermostat related values from parameter file.
+ * \param *fb file buffer containing the config file
+ */
+void config::ParseThermostats(class ConfigFileBuffer * const fb)
+{
+  char * const thermo = new char[12];
+  const int verbose = 0;
+
+  // 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);
+};
+
Index: src/config.hpp
===================================================================
--- src/config.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/config.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -20,28 +20,14 @@
 #include <string>
 
-#include "bondgraph.hpp"
-
 /****************************************** forward declarations *****************************/
 
+class BondGraph;
+class ConfigFileBuffer;
 class molecule;
 class MoleculeListClass;
 class periodentafel;
+class ThermoStatContainer;
 
 /********************************************** declarations *******************************/
-
-class ConfigFileBuffer {
-  public:
-    char **buffer;
-    int *LineMapping;
-    int CurrentLine;
-    int NoLines;
-
-    ConfigFileBuffer();
-    ConfigFileBuffer(const char * const filename);
-    ~ConfigFileBuffer();
-
-    void InitMapping();
-    void MapIonTypesInBuffer(const int NoAtoms);
-};
 
 /** The config file.
@@ -51,4 +37,5 @@
   public:
     class BondGraph *BG;
+    class ThermoStatContainer *Thermostats;
 
     int PsiType;
@@ -70,12 +57,4 @@
     int DoConstrainedMD;
     int MaxOuterStep;
-    int Thermostat;
-    int *ThermostatImplemented;
-    char **ThermostatNames;
-    double TempFrequency;
-    double alpha;
-    double HooverMass;
-    double TargetTemp;
-    int ScaleTempStep;
 
   private:
@@ -152,5 +131,4 @@
   char *GetDefaultPath() const;
   void SetDefaultPath(const char * const path);
-  void InitThermostats();
   void ParseThermostats(class ConfigFileBuffer * const fb);
 };
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/molecule.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -881,5 +881,5 @@
         ElementNo[i] = current++;
     }
-    ActOnAllAtoms( &atom::OutputArrayIndexed, output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
+    ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
     return true;
   }
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/molecule.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -80,8 +80,4 @@
   double *PenaltyConstants;   //!<  penalty constant in front of each term
 };
-
-#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
-enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
-
 
 /** The complete molecule.
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/molecule_dynamics.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -18,4 +18,5 @@
 #include "parser.hpp"
 #include "Plane.hpp"
+#include "ThermoStatContainer.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -643,5 +644,5 @@
 
   // calculate scale configuration
-  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
+  ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
 
   // differentating between the various thermostats
@@ -651,5 +652,5 @@
       break;
      case Woodcock:
-      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
+      if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
         DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
         ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
@@ -684,7 +685,7 @@
       delta_alpha = 0.;
       ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
-      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
-      configuration.alpha += delta_alpha*configuration.Deltat;
-      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
+      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
+      configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
+      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
       // apply updated alpha as additional force
       ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/moleculelist.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -12,4 +12,5 @@
 #include "atom.hpp"
 #include "bond.hpp"
+#include "bondgraph.hpp"
 #include "boundary.hpp"
 #include "config.hpp"
Index: src/unittests/ParserUnitTest.cpp
===================================================================
--- src/unittests/ParserUnitTest.cpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/unittests/ParserUnitTest.cpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -12,6 +12,8 @@
 #include <cppunit/ui/text/TestRunner.h>
 
+#include "Parser/MpqcParser.hpp"
+#include "Parser/PcpParser.hpp"
+#include "Parser/TremoloParser.hpp"
 #include "Parser/XyzParser.hpp"
-#include "Parser/TremoloParser.hpp"
 #include "World.hpp"
 #include "atom.hpp"
@@ -29,10 +31,124 @@
 CPPUNIT_TEST_SUITE_REGISTRATION( ParserUnitTest );
 
+static string waterPcp = "# ParallelCarParinello - main configuration file - created with molecuilder\n\
+\n\
+mainname\tpcp\t# programm name (for runtime files)\n\
+defaultpath\not specified\t# where to put files during runtime\n\
+pseudopotpath\not specified\t# where to find pseudopotentials\n\
+\n\
+ProcPEGamma\t8\t# for parallel computing: share constants\n\
+ProcPEPsi\t1\t# for parallel computing: share wave functions\n\
+DoOutVis\t0\t# Output data for OpenDX\n\
+DoOutMes\t1\t# Output data for measurements\n\
+DoOutOrbitals\t0\t# Output all Orbitals\n\
+DoOutCurr\t0\t# Ouput current density for OpenDx\n\
+DoOutNICS\t0\t# Output Nucleus independent current shieldings\n\
+DoPerturbation\t0\t# Do perturbation calculate and determine susceptibility and shielding\n\
+DoFullCurrent\t0\t# Do full perturbation\n\
+DoConstrainedMD\t0\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD\n\
+Thermostat\tBerendsen\t2.5\t# Which Thermostat and its parameters to use in MD case.\n\
+CommonWannier\t0\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center\n\
+SawtoothStart\t0.01\t# Absolute value for smooth transition at cell border \n\
+VectorPlane\t0\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot\n\
+VectorCut\t0\t# Cut plane axis value\n\
+AddGramSch\t1\t# Additional GramSchmidtOrtogonalization to be safe\n\
+Seed\t1\t# initial value for random seed for Psi coefficients\n\
+\n\
+MaxOuterStep\t0\t# number of MolecularDynamics/Structure optimization steps\n\
+Deltat\t0.01\t# time per MD step\n\
+OutVisStep\t10\t# Output visual data every ...th step\n\
+OutSrcStep\t5\t# Output \"restart\" data every ..th step\n\
+TargetTemp\t0.000950045\t# Target temperature\n\
+MaxPsiStep\t3\t# number of Minimisation steps per state (0 - default)\n\
+EpsWannier\t1e-07\t# tolerance value for spread minimisation of orbitals\n\
+# Values specifying when to stop\n\
+MaxMinStep\t100\t# Maximum number of steps\n\
+RelEpsTotalE\t1e-07\t# relative change in total energy\n\
+RelEpsKineticE\t1e-05\t# relative change in kinetic energy\n\
+MaxMinStopStep\t0\t# check every ..th steps\n\
+MaxMinGapStopStep\t1\t# check every ..th steps\n\
+\n\
+# Values specifying when to stop for INIT, otherwise same as above\n\
+MaxInitMinStep\t100\t# Maximum number of steps\n\
+InitRelEpsTotalE\t1e-05\t# relative change in total energy\n\
+InitRelEpsKineticE\t0.0001\t# relative change in kinetic energy\n\
+InitMaxMinStopStep\t0\t# check every ..th steps\n\
+InitMaxMinGapStopStep\t1\t# check every ..th steps\n\
+\n\
+BoxLength\t# (Length of a unit cell)\n\
+20\n\
+0\t20\n\
+0\t0\t20\n\
+\n\
+ECut\t128\t# energy cutoff for discretization in Hartrees\n\
+MaxLevel\t5\t# number of different levels in the code, >=2\n\
+Level0Factor\t2\t# factor by which node number increases from S to 0 level\n\
+RiemannTensor\t0\t# (Use metric)\n\
+PsiType\t0\t# 0 - doubly occupied, 1 - SpinUp,SpinDown\n\
+MaxPsiDouble\t0\t# here: specifying both maximum number of SpinUp- and -Down-states\n\
+PsiMaxNoUp\t0\t# here: specifying maximum number of SpinUp-states\n\
+PsiMaxNoDown\t0\t# here: specifying maximum number of SpinDown-states\n\
+AddPsis\t0\t# Additional unoccupied Psis for bandgap determination\n\
+\n\
+RCut\t20\t# R-cut for the ewald summation\n\
+StructOpt\t0\t# Do structure optimization beforehand\n\
+IsAngstroem\t1\t# 0 - Bohr, 1 - Angstroem\n\
+RelativeCoord\t0\t# whether ion coordinates are relative (1) or absolute (0)\n\
+MaxTypes\t2\t# maximum number of different ion types\n\
+\n\
+# Ion type data (PP = PseudoPotential, Z = atomic number)\n\
+#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\t# chemical name, symbol\n\
+Ion_Type1\t2\t1\t1.0\t3\t3\t1.008\tHydrogen\tH\n\
+Ion_Type2\t1\t8\t1.0\t3\t3\t15.999\tOxygen\tO\n\
+#Ion_TypeNr._Nr.R[0]\tR[1]\tR[2]\tMoveType (0 MoveIon, 1 FixedIon)\n\
+Ion_Type2_1\t0.000000000\t0.000000000\t0.000000000\t0 # molecule nr 0\n\
+Ion_Type1_1\t0.758602\t0.000000000\t0.504284\t0 # molecule nr 1\n\
+Ion_Type1_2\t0.758602\t0.000000000\t-0.504284\t0 # molecule nr 2\n";
+static string waterMpqc ="% Created by MoleCuilder\n\
+mpqc: (\n\
+\tsavestate = no\n\
+\tdo_gradient = yes\n\
+\tmole<MBPT2>: (\n\
+\t\tmaxiter = 200\n\
+\t\tbasis = $:basis\n\
+\t\tmolecule = $:molecule\n\
+\t\treference<CLHF>: (\n\
+\t\t\tbasis = $:basis\n\
+\t\t\tmolecule = $:molecule\n\
+\t\t)\n\
+\t)\n\
+)\n\
+molecule<Molecule>: (\n\
+\tunit = angstrom\n\
+\t{ atoms geometry } = {\n\
+\t\tO [ 0\t0\t0 ]\n\
+\t\tH [ 0.758602\t0\t0.504284 ]\n\
+\t\tH [ 0.758602\t0\t-0.504284 ]\n\
+\t}\n\
+)\n\
+basis<GaussianBasisSet>: (\n\
+\tname = \"3-21G\"\n\
+\tmolecule = $:molecule\n\
+)\n";
+static string waterXyz = "3\nH2O: water molecule\nO\t0.000000\t0.000000\t0.000000\nH\t0.758602\t0.000000\t0.504284\nH\t0.758602\t0.000000\t-0.504284\n";
+static string Tremolo_Atomdata1 = "# ATOMDATA\tId\tname\tType\tx=3\n";
+static string Tremolo_Atomdata2 = "#\n#ATOMDATA Id name Type x=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
+static string Tremolo_invalidkey = "#\n#ATOMDATA Id name foo Type x=3\n\n\n";
+static string Tremolo_velocity = "#\n#ATOMDATA Id name Type u=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
+static string Tremolo_neighbours = "#\n#ATOMDATA Id Type neighbors=2\n1 H 3 0\n2 H 3 0\n3 O 1 2\n";
+static string Tremolo_improper = "#\n#ATOMDATA Id Type imprData\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
+static string Tremolo_torsion = "#\n#ATOMDATA Id Type torsion\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
+static string Tremolo_full = "# ATOMDATA\tx=3\tu=3\tF\tstress\tId\tneighbors=5\timprData\tGroupMeasureTypeNo\tType\textType\tname\tresName\tchainID\tresSeq\toccupancy\ttempFactor\tsegID\tCharge\tcharge\tGrpTypeNo\ttorsion\n0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t-\t0\tH\t-\t-\t-\t0\t0\t0\t0\t0\t0\t0\t0\t-\t\n";
 
 void ParserUnitTest::setUp() {
   World::getInstance();
+
+  // we need hydrogens and oxygens in the following tests
+  CPPUNIT_ASSERT(World::getInstance().getPeriode()->FindElement(1) != NULL);
+  CPPUNIT_ASSERT(World::getInstance().getPeriode()->FindElement(8) != NULL);
 }
 
 void ParserUnitTest::tearDown() {
+  ChangeTracker::purgeInstance();
   World::purgeInstance();
 }
@@ -43,5 +159,4 @@
   cout << "Testing the XYZ parser." << endl;
   XyzParser* testParser = new XyzParser();
-  string waterXyz = "3\nH2O: water molecule\nO\t0.000000\t0.000000\t0.000000\nH\t0.758602\t0.000000\t0.504284\nH\t0.758602\t0.000000\t-0.504284\n";
   stringstream input;
   input << waterXyz;
@@ -62,18 +177,15 @@
   TremoloParser* testParser = new TremoloParser();
   stringstream input, output;
-  string waterTremolo;
 
   // Atomdata beginning with "# ATOMDATA"
-  waterTremolo = "# ATOMDATA\tId\tname\tType\tx=3\n";
-  input << waterTremolo;
-  testParser->load(&input);
-  testParser->save(&output);
-  CPPUNIT_ASSERT(waterTremolo == output.str());
+  input << Tremolo_Atomdata1;
+  testParser->load(&input);
+  testParser->save(&output);
+  CPPUNIT_ASSERT(Tremolo_Atomdata1 == output.str());
   input.clear();
   output.clear();
 
   // Atomdata beginning with "#ATOMDATA"
-  waterTremolo = "#\n#ATOMDATA Id name Type x=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
-  input << waterTremolo;
+  input << Tremolo_Atomdata2;
   testParser->load(&input);
   testParser->save(&output);
@@ -83,6 +195,5 @@
 
   // Invalid key in Atomdata line
-  waterTremolo = "#\n#ATOMDATA Id name foo Type x=3\n\n\n";
-  input << waterTremolo;
+  input << Tremolo_invalidkey;
   testParser->load(&input);
   //TODO: proove invalidity
@@ -93,9 +204,7 @@
   TremoloParser* testParser = new TremoloParser();
   stringstream input;
-  string waterTremolo;
 
   // One simple data line
-  waterTremolo = "#\n#ATOMDATA Id name Type x=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
-  input << waterTremolo;
+  input << Tremolo_Atomdata2;
   testParser->load(&input);
   CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->x[0] == 3.0);
@@ -106,9 +215,7 @@
   TremoloParser* testParser = new TremoloParser();
   stringstream input;
-  string waterTremolo;
 
   // One simple data line
-  waterTremolo = "#\n#ATOMDATA Id name Type u=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
-  input << waterTremolo;
+  input << Tremolo_velocity;
   testParser->load(&input);
   CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->v[0] == 3.0);
@@ -119,9 +226,7 @@
   TremoloParser* testParser = new TremoloParser();
   stringstream input;
-  string waterTremolo;
 
   // Neighbor data
-  waterTremolo = "#\n#ATOMDATA Id Type neighbors=2\n1 H 3 0\n2 H 3 0\n3 O 1 2\n";
-  input << waterTremolo;
+  input << Tremolo_neighbours;
   testParser->load(&input);
 
@@ -135,9 +240,7 @@
   TremoloParser* testParser = new TremoloParser();
   stringstream input, output;
-  string waterTremolo;
 
   // Neighbor data
-  waterTremolo = "#\n#ATOMDATA Id Type imprData\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
-  input << waterTremolo;
+  input << Tremolo_improper;
   testParser->load(&input);
   testParser->save(&output);
@@ -151,9 +254,7 @@
   TremoloParser* testParser = new TremoloParser();
   stringstream input, output;
-  string waterTremolo;
 
   // Neighbor data
-  waterTremolo = "#\n#ATOMDATA Id Type torsion\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
-  input << waterTremolo;
+  input << Tremolo_torsion;
   testParser->load(&input);
   testParser->save(&output);
@@ -173,6 +274,51 @@
   testParser->setFieldsForSave("x=3 u=3 F stress Id neighbors=5 imprData GroupMeasureTypeNo Type extType name resName chainID resSeq occupancy tempFactor segID Charge charge GrpTypeNo torsion");
   testParser->save(&output);
-  CPPUNIT_ASSERT(output.str() == "# ATOMDATA\tx=3\tu=3\tF\tstress\tId\tneighbors=5\timprData\tGroupMeasureTypeNo\tType\textType\tname\tresName\tchainID\tresSeq\toccupancy\ttempFactor\tsegID\tCharge\tcharge\tGrpTypeNo\ttorsion\n0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t-\t0\tH\t-\t-\t-\t0\t0\t0\t0\t0\t0\t0\t0\t-\t\n");
+  CPPUNIT_ASSERT(output.str() == Tremolo_full);
 
   cout << "testing the tremolo parser is done" << endl;
 }
+
+void ParserUnitTest::readwritePcpTest() {
+  stringstream input(waterPcp);
+  PcpParser* testParser = new PcpParser();
+  testParser->load(&input);
+  input.clear();
+
+  CPPUNIT_ASSERT_EQUAL(3, World::getInstance().numAtoms());
+
+  string newWaterPcp = "";
+  stringstream output;
+  testParser->save(&output);
+
+  input << output;
+  PcpParser* testParser2 = new PcpParser();
+  testParser2->load(&input);
+
+  CPPUNIT_ASSERT_EQUAL(6, World::getInstance().numAtoms());
+
+  CPPUNIT_ASSERT(*testParser == *testParser2);
+}
+
+void ParserUnitTest::writeMpqcTest() {
+  // build up water molecule
+  atom *Walker = NULL;
+  Walker = World::getInstance().createAtom();
+  Walker->type = World::getInstance().getPeriode()->FindElement(8);
+  Walker->x = Vector(0,0,0);
+  Walker = World::getInstance().createAtom();
+  Walker->type = World::getInstance().getPeriode()->FindElement(1);
+  Walker->x = Vector(0.758602,0,0.504284);
+  Walker = World::getInstance().createAtom();
+  Walker->type = World::getInstance().getPeriode()->FindElement(1);
+  Walker->x = Vector(0.758602,0,-0.504284);
+  CPPUNIT_ASSERT_EQUAL(3, World::getInstance().numAtoms());
+
+  // create two stringstreams, one stored, one created
+  stringstream input(waterMpqc);
+  MpqcParser* testParser = new MpqcParser();
+  stringstream output;
+  testParser->save(&output);
+
+  // compare both configs
+  CPPUNIT_ASSERT(input.str() == output.str());
+}
Index: src/unittests/ParserUnitTest.hpp
===================================================================
--- src/unittests/ParserUnitTest.hpp	(revision f34c238ce5cc901d49c2b7d96c15786b1d2d390d)
+++ src/unittests/ParserUnitTest.hpp	(revision d3dff06e8fff78fceae7c93eb68106e91329b667)
@@ -22,4 +22,6 @@
   CPPUNIT_TEST ( readAndWriteTremoloTorsionInformationTest );
   CPPUNIT_TEST ( writeTremoloTest );
+  CPPUNIT_TEST ( readwritePcpTest );
+  CPPUNIT_TEST ( writeMpqcTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -36,4 +38,6 @@
   void readAndWriteTremoloTorsionInformationTest();
   void writeTremoloTest();
+  void readwritePcpTest();
+  void writeMpqcTest();
 };
 
