Changes in src/config.cpp [68f03d:35b698]

File:

• src/config.cpp (modified) (43 diffs)

src/config.cpp

@@ -5 +5 @@
  */
 
+#include "Helpers/MemDebug.hpp"
+
 #include <stdio.h>
 #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"
+#include "info.hpp"
 #include "lists.hpp"
 #include "log.hpp"
@@ -20 +24 @@
 #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 = Malloc<char*>(NoLines, "ConfigFileBuffer::ConfigFileBuffer: **buffer");
-
-  // scan each line and put into buffer
-  int lines=0;
-  int i;
-  do {
-    buffer[lines] = Malloc<char>(MAXSTRINGSIZE, "ConfigFileBuffer::ConfigFileBuffer: *buffer[]");
-    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)
-    Free(&buffer[i]);
-  Free(&buffer);
-  Free(&LineMapping);
-}
-
-
-/** Create trivial mapping.
- */
-void ConfigFileBuffer::InitMapping()
-{
-  LineMapping = Malloc<int>(NoLines, "ConfigFileBuffer::InitMapping: *LineMapping");
-  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 ***************************/
@@ -171 +31 @@
 /** 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),
+    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),
@@ -179 +38 @@
     MaxLevel(5), RiemannTensor(0), LevRFactor(0), RiemannLevel(0), Lev0Factor(2), RTActualUse(0), AddPsis(0), RCut(20.), StructOpt(0), IsAngstroem(1), RelativeCoord(0),
     MaxTypes(0) {
-  mainname = Malloc<char>(MAXSTRINGSIZE,"config constructor: mainname");
-  defaultpath = Malloc<char>(MAXSTRINGSIZE,"config constructor: defaultpath");
-  pseudopotpath = Malloc<char>(MAXSTRINGSIZE,"config constructor: pseudopotpath");
-  databasepath = Malloc<char>(MAXSTRINGSIZE,"config constructor: databasepath");
-  configpath = Malloc<char>(MAXSTRINGSIZE,"config constructor: configpath");
-  configname = Malloc<char>(MAXSTRINGSIZE,"config constructor: configname");
+  mainname = new char[MAXSTRINGSIZE];
+  defaultpath = new char[MAXSTRINGSIZE];
+  pseudopotpath = new char[MAXSTRINGSIZE];
+  databasepath = new char[MAXSTRINGSIZE];
+  configname = new char[MAXSTRINGSIZE];
+  Thermostats = new ThermoStatContainer();
   strcpy(mainname,"pcp");
   strcpy(defaultpath,"not specified");
   strcpy(pseudopotpath,"not specified");
-  configpath[0]='\0';
   configname[0]='\0';
   basis = "3-21G";
-
-  InitThermostats();
 };
 
@@ -199 +55 @@
 config::~config()
 {
-  Free(&mainname);
-  Free(&defaultpath);
-  Free(&pseudopotpath);
-  Free(&databasepath);
-  Free(&configpath);
-  Free(&configname);
-  Free(&ThermostatImplemented);
-  for (int j=0;j<MaxThermostats;j++)
-    Free(&ThermostatNames[j]);
-  Free(&ThermostatNames);
+  delete[](mainname);
+  delete[](defaultpath);
+  delete[](pseudopotpath);
+  delete[](databasepath);
+  delete[](configname);
+  if (Thermostats != NULL)
+    delete(Thermostats);
 
   if (BG != NULL)
     delete(BG);
 };
-
-/** Initialises variables in class config for Thermostats.
- */
-void config::InitThermostats()
-{
-  ThermostatImplemented = Malloc<int>(MaxThermostats, "config constructor: *ThermostatImplemented");
-  ThermostatNames = Malloc<char*>(MaxThermostats, "config constructor: *ThermostatNames");
-  for (int j=0;j<MaxThermostats;j++)
-    ThermostatNames[j] = Malloc<char>(12, "config constructor: ThermostatNames[]");
-
-  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 = Malloc<char>(12, "IonsInitRead: thermo");
-  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;
-  }
-  Free(thermo);
-};
-
 
 /** Displays menu for editing each entry of the config file.
@@ -626 +380 @@
 };
 
-/** Retrieves the path in the given config file name.
- * \param filename config file string
- */
-void config::RetrieveConfigPathAndName(const string filename)
-{
-  char *ptr = NULL;
-  char *buffer = new char[MAXSTRINGSIZE];
-  strncpy(buffer, filename.c_str(), MAXSTRINGSIZE);
-  int last = -1;
-  for(last=MAXSTRINGSIZE;last--;) {
-    if (buffer[last] == '/')
-      break;
-  }
-  if (last == -1) { // no path in front, set to local directory.
-    strcpy(configpath, "./");
-    ptr = buffer;
-  } else {
-    strncpy(configpath, buffer, last+1);
-    ptr = &buffer[last+1];
-    if (last < 254)
-      configpath[last+1]='\0';
-  }
-  strcpy(configname, ptr);
-  DoLog(0) && (Log() << Verbose(0) << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl);
-  delete[](buffer);
-};
-
-/** 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
@@ -861 +573 @@
   file->close();
   delete(file);
-  RetrieveConfigPathAndName(filename);
 
   // ParseParameterFile
-  struct ConfigFileBuffer *FileBuffer = NULL;
-  PrepareFileBuffer(filename,FileBuffer);
+  class ConfigFileBuffer *FileBuffer = new ConfigFileBuffer(filename);
 
   /* Oeffne Hauptparameterdatei */
@@ -874 +584 @@
   int verbose = 0;
   
+  //TODO: This is actually sensible?: if (MaxOuterStep > 0)
   ParseThermostats(FileBuffer);
   
@@ -938 +649 @@
   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))
@@ -1098 +809 @@
     return;
   }
-  RetrieveConfigPathAndName(filename);
   // ParseParameters
 
@@ -1147 +857 @@
   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;
 
@@ -1336 +1046 @@
     *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;
     };
@@ -1369 +1079 @@
     *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;
@@ -1480 +1190 @@
     // 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;
@@ -1522 +1232 @@
     // 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;
@@ -1547 +1257 @@
   int AtomNo = -1;
   int MolNo = 0;
-  atom *Walker = NULL;
   FILE *f = NULL;
 
@@ -1560 +1269 @@
   fprintf(f, "# Created by MoleCuilder\n");
 
-  for (MoleculeList::const_iterator Runner = MolList->ListOfMolecules.begin(); Runner != MolList->ListOfMolecules.end(); Runner++) {
-    Walker = (*Runner)->start;
-    int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
+  for (MoleculeList::const_iterator MolRunner = MolList->ListOfMolecules.begin(); MolRunner != MolList->ListOfMolecules.end(); MolRunner++) {
+    int *elementNo = new int[MAX_ELEMENTS];
+    for (int i=0;i<MAX_ELEMENTS;i++)
+      elementNo[i] = 0;
     AtomNo = 0;
-    while (Walker->next != (*Runner)->end) {
-      Walker = Walker->next;
-      sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
-      elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
+    for (molecule::const_iterator iter = (*MolRunner)->begin(); iter != (*MolRunner)->end(); ++iter) {
+      sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
+      elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
       fprintf(f,
              "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
-             Walker->nr,                /* atom serial number */
+             (*iter)->nr,                /* atom serial number */
              name,         /* atom name */
-             (*Runner)->name,      /* residue name */
+             (*MolRunner)->name,      /* residue name */
              'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
              MolNo,         /* residue sequence number */
-             Walker->node->at(0),                 /* position X in Angstroem */
-             Walker->node->at(1),                 /* position Y in Angstroem */
-             Walker->node->at(2),                 /* position Z in Angstroem */
-             (double)Walker->type->Valence,         /* occupancy */
-             (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
+             (*iter)->node->at(0),                 /* position X in Angstroem */
+             (*iter)->node->at(1),                 /* position Y in Angstroem */
+             (*iter)->node->at(2),                 /* position Z in Angstroem */
+             (double)(*iter)->type->Valence,         /* occupancy */
+             (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
              "0",            /* segment identifier */
-             Walker->type->symbol,    /* element symbol */
+             (*iter)->type->symbol,    /* element symbol */
              "0");           /* charge */
       AtomNo++;
     }
-    Free(&elementNo);
+    delete[](elementNo);
     MolNo++;
   }
@@ -1601 +1310 @@
 {
   int AtomNo = -1;
-  atom *Walker = NULL;
   FILE *f = NULL;
 
-  int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
+  int *elementNo = new int[MAX_ELEMENTS];
+  for (int i=0;i<MAX_ELEMENTS;i++)
+    elementNo[i] = 0;
   char name[MAXSTRINGSIZE];
   strncpy(name, filename, MAXSTRINGSIZE-1);
@@ -1611 +1321 @@
   if (f == NULL) {
     DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
-    Free(&elementNo);
+    delete[](elementNo);
     return false;
   }
   fprintf(f, "# Created by MoleCuilder\n");
 
-  Walker = mol->start;
   AtomNo = 0;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
-    elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
+    elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
     fprintf(f,
            "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
-           Walker->nr,                /* atom serial number */
+           (*iter)->nr,                /* atom serial number */
            name,         /* atom name */
            mol->name,      /* residue name */
            'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
            0,         /* residue sequence number */
-           Walker->node->at(0),                 /* position X in Angstroem */
-           Walker->node->at(1),                 /* position Y in Angstroem */
-           Walker->node->at(2),                 /* position Z in Angstroem */
-           (double)Walker->type->Valence,         /* occupancy */
-           (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
+           (*iter)->node->at(0),                 /* position X in Angstroem */
+           (*iter)->node->at(1),                 /* position Y in Angstroem */
+           (*iter)->node->at(2),                 /* position Z in Angstroem */
+           (double)(*iter)->type->Valence,         /* occupancy */
+           (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
            "0",            /* segment identifier */
-           Walker->type->symbol,    /* element symbol */
+           (*iter)->type->symbol,    /* element symbol */
            "0");           /* charge */
     AtomNo++;
   }
   fclose(f);
-  Free(&elementNo);
+  delete[](elementNo);
 
   return true;
@@ -1653 +1361 @@
 bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
 {
-  atom *Walker = NULL;
   ofstream *output = NULL;
   stringstream * const fname = new stringstream;
@@ -1666 +1373 @@
 
   // scan maximum number of neighbours
-  Walker = mol->start;
   int MaxNeighbours = 0;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    const int count = Walker->ListOfBonds.size();
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    const int count = (*iter)->ListOfBonds.size();
     if (MaxNeighbours < count)
       MaxNeighbours = count;
   }
-  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
-
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    *output << Walker->nr << "\t";
-    *output << Walker->getName() << "\t";
+  *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
+
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    *output << (*iter)->nr << "\t";
+    *output << (*iter)->getName() << "\t";
     *output << mol->name << "\t";
     *output << 0 << "\t";
-    *output << Walker->node->at(0) << "\t" << Walker->node->at(1) << "\t" << Walker->node->at(2) << "\t";
-    *output << static_cast<double>(Walker->type->Valence) << "\t";
-    *output << Walker->type->symbol << "\t";
-    for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
-      *output << (*runner)->GetOtherAtom(Walker)->nr << "\t";
-    for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
+    *output << (*iter)->node->at(0) << "\t" << (*iter)->node->at(1) << "\t" << (*iter)->node->at(2) << "\t";
+    *output << static_cast<double>((*iter)->type->Valence) << "\t";
+    *output << (*iter)->type->symbol << "\t";
+    for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
+      *output << (*runner)->GetOtherAtom(*iter)->nr << "\t";
+    for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
       *output << "-\t";
     *output << endl;
@@ -1708 +1411 @@
 bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
 {
-  atom *Walker = NULL;
+  Info FunctionInfo(__func__);
   ofstream *output = NULL;
   stringstream * const fname = new stringstream;
@@ -1723 +1426 @@
   int MaxNeighbours = 0;
   for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
-    Walker = (*MolWalker)->start;
-    while (Walker->next != (*MolWalker)->end) {
-      Walker = Walker->next;
-      const int count = Walker->ListOfBonds.size();
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      const int count = (*iter)->ListOfBonds.size();
       if (MaxNeighbours < count)
         MaxNeighbours = count;
     }
   }
-  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
+  *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
 
   // create global to local id map
-  int **LocalNotoGlobalNoMap = Calloc<int *>(MolList->ListOfMolecules.size(), "config::SaveTREMOLO - **LocalNotoGlobalNoMap");
+  map<int, int> LocalNotoGlobalNoMap;
   {
-    int MolCounter = 0;
-    int AtomNo = 0;
+    unsigned int MolCounter = 0;
+    int AtomNo = 1;
     for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
-      LocalNotoGlobalNoMap[MolCounter] = Calloc<int>(MolList->CountAllAtoms(), "config::SaveTREMOLO - *LocalNotoGlobalNoMap[]");
-
-      (*MolWalker)->SetIndexedArrayForEachAtomTo( LocalNotoGlobalNoMap[MolCounter], &atom::nr, IncrementalAbsoluteValue, &AtomNo);
-
+      for(molecule::iterator AtomRunner = (*MolWalker)->begin(); AtomRunner != (*MolWalker)->end(); ++AtomRunner) {
+        LocalNotoGlobalNoMap.insert( pair<int,int>((*AtomRunner)->getId(), AtomNo++) );
+      }
       MolCounter++;
     }
+    ASSERT(MolCounter == MolList->ListOfMolecules.size(), "SaveTREMOLO: LocalNotoGlobalNoMap[] has not been correctly initialized for each molecule");
   }
 
@@ -1752 +1453 @@
     int AtomNo = 0;
     for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
-      Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        *output << AtomNo+1 << "\t";
-        *output << Walker->getName() << "\t";
+      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+        *output << LocalNotoGlobalNoMap[ (*iter)->getId() ] << "\t";
+        *output << (*iter)->getName() << "\t";
         *output << (*MolWalker)->name << "\t";
         *output << MolCounter+1 << "\t";
-        *output << Walker->node->at(0) << "\t" << Walker->node->at(1) << "\t" << Walker->node->at(2) << "\t";
-        *output << (double)Walker->type->Valence << "\t";
-        *output << Walker->type->symbol << "\t";
-        for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
-          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ]+1 << "\t";
-        for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
+        *output << (*iter)->node->at(0) << "\t" << (*iter)->node->at(1) << "\t" << (*iter)->node->at(2) << "\t";
+        *output << (double)(*iter)->type->Valence << "\t";
+        *output << (*iter)->type->symbol << "\t";
+        for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
+          *output << LocalNotoGlobalNoMap[ (*runner)->GetOtherAtom((*iter))->getId() ] << "\t";
+        for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
           *output << "-\t";
         *output << endl;
@@ -1778 +1477 @@
   delete(output);
   delete(fname);
-  for(size_t i=0;i<MolList->ListOfMolecules.size(); i++)
-    Free(&LocalNotoGlobalNoMap[i]);
-  Free(&LocalNotoGlobalNoMap);
 
   return true;
@@ -1795 +1491 @@
   char filename[MAXSTRINGSIZE];
   ofstream output;
-  molecule *mol = World::getInstance().createMolecule();
-  mol->SetNameFromFilename(ConfigFileName);
-
-  if (!strcmp(configpath, GetDefaultPath())) {
-    eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
-  }
-
+  molecule *mol = NULL;
 
   // first save as PDB data
@@ -1808 +1498 @@
   if (output == NULL)
     strcpy(filename,"main_pcp_linux");
-  Log() << Verbose(0) << "Saving as pdb input ";
+  Log() << Verbose(0) << "Saving as pdb input ... " << endl;
   if (SavePDB(filename, molecules))
-    Log() << Verbose(0) << "done." << endl;
+    Log() << Verbose(0) << "\t... done." << endl;
   else
-    Log() << Verbose(0) << "failed." << endl;
+    Log() << Verbose(0) << "\t... failed." << endl;
 
   // then save as tremolo data file
@@ -1819 +1509 @@
   if (output == NULL)
     strcpy(filename,"main_pcp_linux");
-  Log() << Verbose(0) << "Saving as tremolo data input ";
+  Log() << Verbose(0) << "Saving as tremolo data input ... " << endl;
   if (SaveTREMOLO(filename, molecules))
-    Log() << Verbose(0) << "done." << endl;
+    Log() << Verbose(0) << "\t... done." << endl;
   else
-    Log() << Verbose(0) << "failed." << endl;
+    Log() << Verbose(0) << "\t... failed." << endl;
 
   // translate each to its center and merge all molecules in MoleculeListClass into this molecule
   int N = molecules->ListOfMolecules.size();
-  int *src = new int[N];
-  N=0;
-  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
-    src[N++] = (*ListRunner)->IndexNr;
-    (*ListRunner)->Translate(&(*ListRunner)->Center);
-  }
-  molecules->SimpleMultiAdd(mol, src, N);
-  delete[](src);
-
-  // ... and translate back
-  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
-    (*ListRunner)->Center.Scale(-1.);
-    (*ListRunner)->Translate(&(*ListRunner)->Center);
-    (*ListRunner)->Center.Scale(-1.);
+  if (N != 1) { // don't do anything in case of only one molecule (shifts mol ids otherwise)
+    int *src = new int[N];
+    N=0;
+    for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      src[N++] = (*ListRunner)->IndexNr;
+      (*ListRunner)->Translate(&(*ListRunner)->Center);
+    }
+    mol = World::getInstance().createMolecule();
+    mol->SetNameFromFilename(ConfigFileName);
+    molecules->SimpleMultiMerge(mol, src, N);
+    mol->doCountAtoms();
+    mol->CountElements();
+    //mol->CalculateOrbitals(*this);
+    delete[](src);
+  } else {
+    if (!molecules->ListOfMolecules.empty()) {
+      mol = *(molecules->ListOfMolecules.begin());
+      mol->doCountAtoms();
+      //mol->CalculateOrbitals(*this);
+    } else {
+      DoeLog(1) && (eLog() << Verbose(1) << "There are no molecules to save!" << endl);
+    }
   }
 
   Log() << Verbose(0) << "Storing configuration ... " << endl;
   // get correct valence orbitals
-  mol->CalculateOrbitals(*this);
-  InitMaxMinStopStep = MaxMinStopStep = MaxPsiDouble;
   if (ConfigFileName != NULL) { // test the file name
     strcpy(filename, ConfigFileName);
@@ -1859 +1555 @@
   output.close();
   output.clear();
-  Log() << Verbose(0) << "Saving of config file ";
+  Log() << Verbose(0) << "Saving of config file ... " << endl;
   if (Save(filename, periode, mol))
-    Log() << Verbose(0) << "successful." << endl;
+    Log() << Verbose(0) << "\t... successful." << endl;
   else
-    Log() << Verbose(0) << "failed." << endl;
+    Log() << Verbose(0) << "\t... failed." << endl;
 
   // and save to xyz file
@@ -1876 +1572 @@
     output.open(filename, ios::trunc);
   }
-  Log() << Verbose(0) << "Saving of XYZ file ";
+  Log() << Verbose(0) << "Saving of XYZ file ... " << endl;
   if (mol->MDSteps <= 1) {
     if (mol->OutputXYZ(&output))
-      Log() << Verbose(0) << "successful." << endl;
+      Log() << Verbose(0) << "\t... successful." << endl;
     else
-      Log() << Verbose(0) << "failed." << endl;
+      Log() << Verbose(0) << "\t... failed." << endl;
   } else {
     if (mol->OutputTrajectoriesXYZ(&output))
-      Log() << Verbose(0) << "successful." << endl;
+      Log() << Verbose(0) << "\t... successful." << endl;
     else
-      Log() << Verbose(0) << "failed." << endl;
+      Log() << Verbose(0) << "\t... failed." << endl;
   }
   output.close();
@@ -1896 +1592 @@
   if (output == NULL)
     strcpy(filename,"main_pcp_linux");
-  Log() << Verbose(0) << "Saving as mpqc input ";
+  Log() << Verbose(0) << "Saving as mpqc input .. " << endl;
   if (SaveMPQC(filename, mol))
-    Log() << Verbose(0) << "done." << endl;
+    Log() << Verbose(0) << "\t... done." << endl;
   else
-    Log() << Verbose(0) << "failed." << endl;
-
-  if (!strcmp(configpath, GetDefaultPath())) {
-    eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
-  }
-
-  World::getInstance().destroyMolecule(mol);
+    Log() << Verbose(0) << "\t... failed." << endl;
+
+  // don't destroy molecule as it contains all our atoms
+  //World::getInstance().destroyMolecule(mol);
 };
 
@@ -1938 +1631 @@
   char *dummy1 = NULL;
   char *dummy = NULL;
-  char * const free_dummy = Malloc<char>(256, "config::ParseForParameter: *free_dummy");    // pointers in the line that is read in per step
+  char free_dummy[MAXSTRINGSIZE];    // pointers in the line that is read in per step
   dummy1 = free_dummy;
 
@@ -1954 +1647 @@
       if (file->eof()) {
         if ((critical) && (found == 0)) {
-          Free(free_dummy);
           //Error(InitReading, name);
           fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
@@ -1962 +1654 @@
           file->clear();
           file->seekg(file_position, ios::beg);  // rewind to start position
-          Free(free_dummy);
           return 0;
         }
@@ -1993 +1684 @@
         dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
         //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
-        //Free((void **)&free_dummy);
         //Error(FileOpenParams, NULL);
       } else {
@@ -2014 +1704 @@
               if (file->eof()) {
                 if ((critical) && (found == 0)) {
-                  Free(free_dummy);
                   //Error(InitReading, name);
                   fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
@@ -2022 +1711 @@
                   file->clear();
                   file->seekg(file_position, ios::beg);  // rewind to start position
-                  Free(free_dummy);
                   return 0;
                 }
@@ -2063 +1751 @@
                   if (critical) {
                     if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
-                    Free(free_dummy);
                     //return 0;
                     exit(255);
@@ -2071 +1758 @@
                     file->clear();
                     file->seekg(file_position, ios::beg);  // rewind to start position
-                    Free(free_dummy);
                     return 0;
                   }
@@ -2084 +1770 @@
                   file->seekg(file_position, ios::beg);  // rewind to start position
                 }
-                Free(free_dummy);
                 return 0;
               }
@@ -2140 +1825 @@
   if ((type >= row_int) && (verbose))
     fprintf(stderr,"\n");
-  Free(free_dummy);
   if (!sequential) {
     file->clear();
@@ -2221 +1905 @@
         dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
         //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
-        //Free(&free_dummy);
         //Error(FileOpenParams, NULL);
       } else {
@@ -2361 +2044 @@
   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);
+};
+