Index: src/Parser/PcpParser.cpp
===================================================================
--- src/Parser/PcpParser.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Parser/PcpParser.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -402,26 +402,6 @@
     *file << "DoConstrainedMD\t" << Switches.DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
     ASSERT(Thermostats != NULL, "PcpParser::save() - Thermostats not initialized!");
-    ASSERT(Thermostats->ThermostatNames != NULL, "PcpParser::save() - Thermostats not initialized!");
-    *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 << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
+    *file << Thermostats->activeThermostat->writeParams();
     *file << "\t# Which Thermostat and its parameters to use in MD case." << endl;
     *file << "CommonWannier\t" << LocalizedOrbitals.CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
@@ -588,65 +568,9 @@
   // 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;
-    }
+    Thermostats->makeActive(thermo,fb);
   } else {
     if ((Thermostats->TargetTemp != 0))
       DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
-    Thermostats->Thermostat = None;
+    Thermostats->chooseNone();
   }
   delete[](thermo);
Index: src/ThermoStatContainer.cpp
===================================================================
--- src/ThermoStatContainer.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/ThermoStatContainer.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -21,34 +21,70 @@
 #include "Helpers/Verbose.hpp"
 
+#include <Thermostats/Berendsen.hpp>
+#include <Thermostats/GaussianThermostat.hpp>
+#include <Thermostats/Langevin.hpp>
+#include <Thermostats/NoseHoover.hpp>
+#include <Thermostats/NoThermostat.hpp>
+#include <Thermostats/Woodcock.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)
+  activeThermostat(0),
+  TargetTemp(0.00095004455)
 {
-  ThermostatImplemented = new int[MaxThermostats];
-  ThermostatNames = new char *[MaxThermostats];
-  for (int j=0;j<MaxThermostats;j++)
-    ThermostatNames[j] = new char[12];
+  ThermostatTraits<Thermostat> *BerendsenTrait = new ThermostatTraits<Berendsen>();
+  availThermostats[BerendsenTrait->getName()] = BerendsenTrait;
+  ThermostatTraits<Thermostat> *GaussianTrait = new ThermostatTraits<GaussianThermostat>();
+  availThermostats[GaussianTrait->getName()] = GaussianTrait;
+  ThermostatTraits<Thermostat> *LangevinTrait = new ThermostatTraits<Langevin>();
+  availThermostats[LangevinTrait->getName()] = LangevinTrait;
+  ThermostatTraits<Thermostat> *NoseHooverTrait = new ThermostatTraits<NoseHoover>();
+  availThermostats[NoseHooverTrait->getName()] = NoseHooverTrait;
+  ThermostatTraits<Thermostat> *NoThermostatTrait = new ThermostatTraits<NoThermostat>();
+  availThermostats[NoThermostatTrait->getName()] = NoThermostatTrait;
+  ThermostatTraits<Thermostat> *WoodcockTrait = new ThermostatTraits<Woodcock>();
+  availThermostats[WoodcockTrait->getName()] = WoodcockTrait;
+  for(traitsMap::iterator iter = availThermostats.begin();iter!=availThermostats.end();++iter){
+    cout << iter->first << endl;
+  }
+  ASSERT(availThermostats.size()==6,"Not all implemented thermostats referenced!\nDid you check the names in the traits?");
+  activeThermostat=new Berendsen();
+}
 
-  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;
+ThermostatTraits<Thermostat> *ThermoStatContainer::getTraitByName(const std::string name){
+  if(!availThermostats.count(name))
+    return 0;
+  return availThermostats[name];
+}
+
+Thermostat *ThermoStatContainer::makeByName(const std::string name,class ConfigFileBuffer * const fb){
+  ThermostatTraits<Thermostat>* trait = getTraitByName(name);
+  if(trait){
+    Thermostat *res = trait->make(fb);
+    res->addToContainer(this);
+    return res;
+  }
+  else{
+    return 0;
+  }
+}
+
+void ThermoStatContainer::makeActive(const std::string name,class ConfigFileBuffer * const fb){
+  Thermostat* newThermostat = makeByName(name,fb);
+  if(newThermostat){
+    if(activeThermostat){
+      delete activeThermostat;
+    }
+    activeThermostat = newThermostat;
+  }
+}
+
+void ThermoStatContainer::chooseNone(){
+  if(activeThermostat){
+    delete activeThermostat;
+  }
+  activeThermostat = new NoThermostat();
 }
 
@@ -58,8 +94,10 @@
 ThermoStatContainer::~ThermoStatContainer()
 {
-  delete[](ThermostatImplemented);
-  for (int j=0;j<MaxThermostats;j++)
-    delete[](ThermostatNames[j]);
-  delete[](ThermostatNames);
+  if(activeThermostat){
+    delete activeThermostat;
+  }
+  for(traitsMap::iterator iter= availThermostats.begin();iter!=availThermostats.end();++iter){
+    delete (iter->second);
+  }
 }
 
Index: src/ThermoStatContainer.hpp
===================================================================
--- src/ThermoStatContainer.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/ThermoStatContainer.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -14,8 +14,7 @@
 #endif
 
+#include "Thermostats/Thermostat.hpp"
+
 class ConfigFileBuffer;
-
-enum { 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 {
@@ -23,13 +22,13 @@
   ThermoStatContainer();
   ~ThermoStatContainer();
+  ThermostatTraits<Thermostat> *getTraitByName(const std::string name);
+  Thermostat *makeByName(const std::string name,class ConfigFileBuffer * const fb);
+  void makeActive(const std::string name,class ConfigFileBuffer * const fb);
+  void chooseNone();
 
-  int Thermostat;
-  int *ThermostatImplemented;
-  char **ThermostatNames;
-  double TempFrequency;
-  double alpha;
-  double HooverMass;
+  Thermostat* activeThermostat;
+  typedef std::map<std::string,ThermostatTraits<Thermostat>* > traitsMap;
+  traitsMap availThermostats;
   double TargetTemp;
-  int ScaleTempStep;
 };
 
Index: src/Thermostats/Berendsen.cpp
===================================================================
--- src/Thermostats/Berendsen.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Berendsen.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -25,13 +25,15 @@
 {}
 
-ThermostatTraits<class Berendsen>::ThermostatTraits() :
-  name("Berendsen")
-{}
+const char *ThermostatTraits<Berendsen>::name = "Berendsen";
 
-Thermostat *ThermostatTraits<class Berendsen>::make(class ConfigFileBuffer * const fb){
+std::string ThermostatTraits<Berendsen>::getName(){
+  return ThermostatTraits<Berendsen>::name;
+}
+
+Thermostat *ThermostatTraits<Berendsen>::make(class ConfigFileBuffer * const fb){
   double TempFrequency;
   const int verbose = 0;
   ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
-  return new class Berendsen(TempFrequency);
+  return new Berendsen(TempFrequency);
 }
 
@@ -51,5 +53,5 @@
 double Berendsen::doScaleAtoms(unsigned int step,double ActualTemp,ForwardIterator begin, ForwardIterator end){
   DoLog(2) && (Log() << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl);
-  double ekin;
+  double ekin=0;
   double ScaleTempFactor = getContainer().TargetTemp/ActualTemp;
   for(ForwardIterator iter=begin;iter!=end;++iter){
@@ -64,5 +66,5 @@
 
 std::string Berendsen::name(){
-  return ThermostatTraits<Berendsen>().name;
+  return ThermostatTraits<Berendsen>::name;
 }
 
Index: src/Thermostats/Berendsen.hpp
===================================================================
--- src/Thermostats/Berendsen.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Berendsen.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -36,7 +36,7 @@
 struct ThermostatTraits<class Berendsen> : public ThermostatTraits<Thermostat>
 {
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb);
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/Thermostats/GaussianThermostat.cpp
===================================================================
--- src/Thermostats/GaussianThermostat.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/GaussianThermostat.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -30,7 +30,9 @@
 {}
 
-ThermostatTraits<GaussianThermostat>::ThermostatTraits() :
-  name("Gaussian")
-{}
+const char *ThermostatTraits<GaussianThermostat>::name = "Gaussian";
+
+std::string ThermostatTraits<GaussianThermostat>::getName(){
+  return ThermostatTraits<GaussianThermostat>::name;
+}
 
 Thermostat *ThermostatTraits<GaussianThermostat>::make(class ConfigFileBuffer * const fb){
@@ -93,5 +95,5 @@
 
 std::string GaussianThermostat::name(){
-  return ThermostatTraits<GaussianThermostat>().name;
+  return ThermostatTraits<GaussianThermostat>::name;
 }
 
Index: src/Thermostats/GaussianThermostat.hpp
===================================================================
--- src/Thermostats/GaussianThermostat.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/GaussianThermostat.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -41,7 +41,7 @@
 struct ThermostatTraits<GaussianThermostat> : public ThermostatTraits<Thermostat>
 {
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb);
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/Thermostats/Langevin.cpp
===================================================================
--- src/Thermostats/Langevin.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Langevin.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -32,9 +32,11 @@
 }
 
-ThermostatTraits<class Langevin>::ThermostatTraits() :
-  name("Langevin")
-{}
+const char *ThermostatTraits<Langevin>::name = "Langevin";
 
-Thermostat *ThermostatTraits<class Langevin>::make(class ConfigFileBuffer * const fb){
+std::string ThermostatTraits<Langevin>::getName(){
+  return ThermostatTraits<Langevin>::name;
+}
+
+Thermostat *ThermostatTraits<Langevin>::make(class ConfigFileBuffer * const fb){
   double TempFrequency;
   double alpha;
@@ -46,5 +48,5 @@
     alpha = 1.;
   }
-  return new class Langevin(TempFrequency,alpha);
+  return new Langevin(TempFrequency,alpha);
 }
 
@@ -85,5 +87,5 @@
 
 std::string Langevin::name(){
-  return ThermostatTraits<Langevin>().name;
+  return ThermostatTraits<Langevin>::name;
 }
 
Index: src/Thermostats/Langevin.hpp
===================================================================
--- src/Thermostats/Langevin.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Langevin.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -39,7 +39,7 @@
 struct ThermostatTraits<class Langevin> : public ThermostatTraits<Thermostat>
 {
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb);
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/Thermostats/NoThermostat.cpp
===================================================================
--- src/Thermostats/NoThermostat.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/NoThermostat.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -14,10 +14,12 @@
 {}
 
-ThermostatTraits<class NoThermostat>::ThermostatTraits() :
-  name("None")
-{}
+const char *ThermostatTraits<NoThermostat>::name = "None";
 
-Thermostat *ThermostatTraits<class NoThermostat>::make(class ConfigFileBuffer * const fb){
-  return new class NoThermostat();
+std::string ThermostatTraits<NoThermostat>::getName(){
+  return ThermostatTraits<NoThermostat>::name;
+}
+
+Thermostat *ThermostatTraits<NoThermostat>::make(class ConfigFileBuffer * const fb){
+  return new NoThermostat();
 }
 
@@ -35,5 +37,5 @@
 
 std::string NoThermostat::name(){
-  return ThermostatTraits<NoThermostat>().name;
+  return ThermostatTraits<NoThermostat>::name;
 }
 
Index: src/Thermostats/NoThermostat.hpp
===================================================================
--- src/Thermostats/NoThermostat.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/NoThermostat.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -28,7 +28,7 @@
 struct ThermostatTraits<class NoThermostat>: public ThermostatTraits<Thermostat>
 {
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb);
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/Thermostats/NoseHoover.cpp
===================================================================
--- src/Thermostats/NoseHoover.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/NoseHoover.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -27,13 +27,15 @@
 {}
 
-ThermostatTraits<class NoseHoover>::ThermostatTraits() :
-  name("NoseHoover")
-{}
+const char *ThermostatTraits<NoseHoover>::name = "NoseHoover";
 
-Thermostat *ThermostatTraits<class NoseHoover>::make(class ConfigFileBuffer * const fb){
+std::string ThermostatTraits<NoseHoover>::getName(){
+  return ThermostatTraits<NoseHoover>::name;
+}
+
+Thermostat *ThermostatTraits<NoseHoover>::make(class ConfigFileBuffer * const fb){
   double HooverMass;
   const int verbose = 0;
   ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
-  return new class NoseHoover(HooverMass);
+  return new NoseHoover(HooverMass);
 }
 
@@ -82,5 +84,5 @@
 
 std::string NoseHoover::name(){
-  return ThermostatTraits<NoseHoover>().name;
+  return ThermostatTraits<NoseHoover>::name;
 }
 
Index: src/Thermostats/NoseHoover.hpp
===================================================================
--- src/Thermostats/NoseHoover.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/NoseHoover.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -40,7 +40,7 @@
 struct ThermostatTraits<class NoseHoover> : public ThermostatTraits<Thermostat>
 {
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb);
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/Thermostats/Thermostat.cpp
===================================================================
--- src/Thermostats/Thermostat.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Thermostat.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -17,10 +17,12 @@
 {}
 
-ThermostatTraits<Thermostat>::ThermostatTraits() :
-  name("Base")
-{}
+const char *ThermostatTraits<Thermostat>::name = "Base";
 
-void Thermostat::addToContainer(ThermoStatContainer &_container){
-  container = &_container;
+std::string ThermostatTraits<Thermostat>::getName(){
+  return ThermostatTraits<Thermostat>::name;
+}
+
+void Thermostat::addToContainer(ThermoStatContainer *_container){
+  container = _container;
 }
 
Index: src/Thermostats/Thermostat.hpp
===================================================================
--- src/Thermostats/Thermostat.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Thermostat.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -33,6 +33,6 @@
 
   virtual std::string writeParams()=0;
+  void addToContainer(ThermoStatContainer *);
 protected:
-  void addToContainer(ThermoStatContainer &);
   ThermoStatContainer &getContainer();
 private:
@@ -45,7 +45,7 @@
 template <>
 struct ThermostatTraits<Thermostat>{
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb)=0;
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/Thermostats/Woodcock.cpp
===================================================================
--- src/Thermostats/Woodcock.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Woodcock.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -27,13 +27,15 @@
 {}
 
-ThermostatTraits<class Woodcock>::ThermostatTraits() :
-  name("Woodcock")
-{}
+const char *ThermostatTraits<Woodcock>::name = "Woodcock";
 
-Thermostat *ThermostatTraits<class Woodcock>::make(class ConfigFileBuffer * const fb){
+std::string ThermostatTraits<Woodcock>::getName(){
+  return ThermostatTraits<Woodcock>::name;
+}
+
+Thermostat *ThermostatTraits<Woodcock>::make(class ConfigFileBuffer * const fb){
   int ScaleTempStep;
   const int verbose = 0;
   ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
-  return new class Woodcock(ScaleTempStep);
+  return new Woodcock(ScaleTempStep);
 }
 
@@ -69,5 +71,5 @@
 
 std::string Woodcock::name(){
-  return ThermostatTraits<Woodcock>().name;
+  return ThermostatTraits<Woodcock>::name;
 }
 
Index: src/Thermostats/Woodcock.hpp
===================================================================
--- src/Thermostats/Woodcock.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/Thermostats/Woodcock.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -36,7 +36,7 @@
 struct ThermostatTraits<class Woodcock> : public ThermostatTraits<Thermostat>
 {
-  ThermostatTraits();
   virtual Thermostat *make(class ConfigFileBuffer * const fb);
-  const char* name;
+  virtual std::string getName();
+  static const char* name;
 };
 
Index: src/atom_trajectoryparticle.cpp
===================================================================
--- src/atom_trajectoryparticle.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/atom_trajectoryparticle.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -144,124 +144,4 @@
 };
 
-/** Scales velocity of atom according to Woodcock thermostat.
- * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
- * \param Step MD step to scale
- * \param *ekin sum of kinetic energy
- */
-void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
-{
-  Vector &U = Trajectory.U.at(Step);
-  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      U[d] *= ScaleTempFactor;
-      *ekin += 0.5*getType()->mass * U[d]*U[d];
-    }
-};
-
-/** Scales velocity of atom according to Gaussian thermostat.
- * \param Step MD step to scale
- * \param *G
- * \param *E
- */
-void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
-{
-  Vector &U = Trajectory.U.at(Step);
-  Vector &F = Trajectory.F.at(Step);
-  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      *G += U[d] * F[d];
-      *E += U[d]*U[d]*getType()->mass;
-    }
-};
-
-/** Determines scale factors according to Gaussian thermostat.
- * \param Step MD step to scale
- * \param GE G over E ratio
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and TargetTemp
- */
-void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
-{
-  Vector &U = Trajectory.U.at(Step);
-  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      U[d] += configuration->Deltat/getType()->mass * ( (G_over_E) * (U[d]*getType()->mass) );
-      *ekin += getType()->mass * U[d]*U[d];
-    }
-};
-
-/** Scales velocity of atom according to Langevin thermostat.
- * \param Step MD step to scale
- * \param *r random number generator
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and TargetTemp
- */
-void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
-{
-  double sigma  = sqrt(configuration->Thermostats->TargetTemp/getType()->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->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
-      for (int d=0; d<NDIM; d++) {
-        U[d] = gsl_ran_gaussian (r, sigma);
-      }
-      DoLog(2) && (Log() << Verbose(2) << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl);
-    }
-    for (int d=0; d<NDIM; d++)
-      *ekin += 0.5*getType()->mass * U[d]*U[d];
-  }
-};
-
-/** Scales velocity of atom according to Berendsen thermostat.
- * \param Step MD step to scale
- * \param ScaleTempFactor factor to scale energy (not velocity!) with
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and Deltat
- */
-void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
-{
-  Vector &U = Trajectory.U.at(Step);
-  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->Thermostats->TempFrequency)*(ScaleTempFactor-1));
-      *ekin += 0.5*getType()->mass * U[d]*U[d];
-    }
-  }
-};
-
-/** Initializes current run of NoseHoover thermostat.
- * \param Step MD step to scale
- * \param *delta_alpha additional sum of kinetic energy on return
- */
-void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
-{
-  Vector &U = Trajectory.U.at(Step);
-  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      *delta_alpha += U[d]*U[d]*getType()->mass;
-    }
-  }
-};
-
-/** Initializes current run of NoseHoover thermostat.
- * \param Step MD step to scale
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and Deltat
- */
-void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
-{
-  Vector &U = Trajectory.U.at(Step);
-  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-        U[d] += configuration->Deltat/getType()->mass * (configuration->Thermostats->alpha * (U[d] * getType()->mass));
-        *ekin += (0.5*getType()->mass) * U[d]*U[d];
-      }
-  }
-};
-
-
 std::ostream & TrajectoryParticle::operator << (std::ostream &ost) const
 {
Index: src/atom_trajectoryparticle.hpp
===================================================================
--- src/atom_trajectoryparticle.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/atom_trajectoryparticle.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -54,13 +54,4 @@
   double getKineticEnergy(unsigned int step) const;
 
-  // thermostats
-  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
-  void Thermostat_Gaussian_init(int Step, double *G, double *E);
-  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
-  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
-  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
-  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
-  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
-
   std::ostream & operator << (std::ostream &ost) const;
 
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/config.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -622,5 +622,4 @@
 };
 
-
 /** Stores all elements of config structure from which they can be re-read.
  * \param *filename name of file
@@ -650,25 +649,6 @@
     *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" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
-    switch(Thermostats->Thermostat) {
-      default:
-      case None:
-        break;
-      case Woodcock:
-        *output << Thermostats->ScaleTempStep;
-        break;
-      case Gaussian:
-        *output << Thermostats->ScaleTempStep;
-        break;
-      case Langevin:
-        *output << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
-        break;
-      case Berendsen:
-        *output << Thermostats->TempFrequency;
-        break;
-      case NoseHoover:
-        *output << Thermostats->HooverMass;
-        break;
-    };
+    *output << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
+    *output << Thermostats->activeThermostat->writeParams();
     *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
     *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
@@ -1655,65 +1635,9 @@
   // 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;
-    }
+    Thermostats->makeActive(thermo,fb);
   } else {
     if ((Thermostats->TargetTemp != 0))
       DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
-    Thermostats->Thermostat = None;
+    Thermostats->chooseNone();
   }
   delete[](thermo);
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/molecule.hpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -266,5 +266,4 @@
   void ScanForPeriodicCorrection();
   bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
-  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
   double VolumeOfConvexEnvelope(bool IsAngstroem);
 
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision 8d1dd4fcc7d507ccaf7345578b875e0a40064ac1)
+++ src/molecule_dynamics.cpp	(revision 14c57a9b8ce2c94b3ca4d4d59b36bd5aa57cabf1)
@@ -24,4 +24,5 @@
 #include "LinearAlgebra/Plane.hpp"
 #include "ThermoStatContainer.hpp"
+#include "Thermostats/Berendsen.hpp"
 
 #include <gsl/gsl_matrix.h>
@@ -615,5 +616,8 @@
   ActualTemp = 0.;
   ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
-  Thermostats(configuration, ActualTemp, Berendsen);
+  Berendsen berendsen = Berendsen();
+  berendsen.addToContainer(configuration.Thermostats);
+  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
+  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
   MDSteps++;
 
@@ -621,82 +625,2 @@
   return true;
 };
-
-/** Implementation of various thermostats.
- * All these thermostats apply an additional force which has the following forms:
- * -# Woodcock
- *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
- * -# Gaussian
- *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
- * -# Langevin
- *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$
- * -# Berendsen
- *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
- * -# Nose-Hoover
- *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
- * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
- * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
- * belatedly and the constraint force set.
- * \param *P Problem at hand
- * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
- * \sa InitThermostat()
- */
-void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
-{
-  double ekin = 0.;
-  double E = 0., G = 0.;
-  double delta_alpha = 0.;
-  double ScaleTempFactor;
-  gsl_rng * r;
-  const gsl_rng_type * T;
-
-  // calculate scale configuration
-  ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
-
-  // differentating between the various thermostats
-  switch(Thermostat) {
-     case None:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying no thermostat..." << endl);
-      break;
-     case Woodcock:
-      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 );
-      }
-      break;
-     case Gaussian:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Gaussian thermostat..." << endl);
-      ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
-
-      DoLog(1) && (Log() << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl);
-      ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
-
-      break;
-     case Langevin:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Langevin thermostat..." << endl);
-      // init random number generator
-      gsl_rng_env_setup();
-      T = gsl_rng_default;
-      r = gsl_rng_alloc (T);
-      // Go through each ion
-      ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
-      break;
-
-     case Berendsen:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl);
-      ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
-      break;
-
-     case NoseHoover:
-      DoLog(2) && (Log() << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl);
-      // dynamically evolve alpha (the additional degree of freedom)
-      delta_alpha = 0.;
-      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
-      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 );
-      break;
-  }
-  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
-};
