Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 16a52b3c50fff2ccc678808e2588e1174646b8ec)
+++ src/builder.cpp	(revision 62f793aff2f24a09fa44578b18eca5a99c739605)
@@ -1004,5 +1004,5 @@
               SaveFlag = true;
               cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
-              if (!mol->VerletForceIntegration((ofstream *)&cout, argv[argptr], configuration.Deltat, configuration.GetIsAngstroem(), configuration.DoConstrainedMD))
+              if (!mol->VerletForceIntegration((ofstream *)&cout, argv[argptr], configuration))
                 cout << Verbose(2) << "File not found." << endl;
               else {
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision 16a52b3c50fff2ccc678808e2588e1174646b8ec)
+++ src/config.cpp	(revision 62f793aff2f24a09fa44578b18eca5a99c739605)
@@ -18,4 +18,12 @@
   configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
   configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
+  ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
+  ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
+  for (int j=0;j<MaxThermostats;j++)
+    ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
+  Thermostat = 4;
+  alpha = 0.;
+  ScaleTempStep = 25;
+  TempFrequency = 2.5;
   strcpy(mainname,"pcp");
   strcpy(defaultpath,"not specified");
@@ -23,5 +31,18 @@
   configpath[0]='\0';
   configname[0]='\0';
-  
+
+  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;
+
   FastParsing = false;
   ProcPEGamma=8;
@@ -96,4 +117,80 @@
   Free((void **)&configname, "config::~config: *configname");
 };
+
+/** Readin of Thermostat related values from parameter file.
+ * \param *source parameter file
+ */
+void config::InitThermostats(ifstream *source)
+{
+  char *thermo = MallocString(12, "IonsInitRead: thermo");
+  int verbose = 0;
+
+  // read desired Thermostat from file along with needed additional parameters
+  if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
+    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
+      if (ThermostatImplemented[0] == 1) {
+        Thermostat = None;
+      } else {
+        cout << 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,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
+      } else {
+        cout << 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,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
+      } else {
+        cout << 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,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
+        if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
+          cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
+        } else {
+          alpha = 1.;
+        }
+      } else {
+        cout << 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,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
+      } else {
+        cout << 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,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
+        alpha = 0.;
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else {
+      cout << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
+      Thermostat = None;
+    }
+  } else {
+    if ((MaxOuterStep > 0) && (TargetTemp != 0)) 
+      cout << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
+    Thermostat = None;
+  }
+  Free((void **)&thermo, "InitThermostats: thermo");
+};
+
 
 /** Displays menu for editing each entry of the config file.
@@ -465,4 +562,6 @@
   double value[3];
   
+  InitThermostats(file);
+  
   /* Namen einlesen */
 
@@ -1001,4 +1100,26 @@
     *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) {
+      default:
+      case None:
+        break;
+      case Woodcock:
+        *output << ScaleTempStep;
+        break;
+      case Gaussian:
+        *output << ScaleTempStep;
+        break;
+      case Langevin:
+        *output << TempFrequency << "\t" << alpha;
+        break;
+      case Berendsen:
+        *output << TempFrequency;
+        break;
+      case NoseHoover:
+        *output << HooverMass;
+        break;
+    };
+    *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;
     *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 16a52b3c50fff2ccc678808e2588e1174646b8ec)
+++ src/molecules.cpp	(revision 62f793aff2f24a09fa44578b18eca5a99c739605)
@@ -1017,5 +1017,5 @@
  * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
  * \param *out output stream for debugging
- * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$V^{con}(x)\f$)
+ * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
  * \param startstep start configuration (MDStep in molecule::trajectories)
  * \param endstep end configuration (MDStep in molecule::trajectories)
@@ -1479,4 +1479,5 @@
  * \param *out output stream for debugging
  * \param *file filename
+ * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
  * \param delta_t time step width in atomic units
  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
@@ -1485,7 +1486,6 @@
  * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
  */
-bool molecule::VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained)
-{
-  element *runner = elemente->start;
+bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
+{
   atom *walker = NULL;
   int AtomNo;
@@ -1493,5 +1493,5 @@
   string token;
   stringstream item;
-  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy;
+  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp = 0.;
   ForceMatrix Force;
 
@@ -1523,82 +1523,242 @@
       }
     // solve a constrained potential if we are meant to
-    if (DoConstrained) {
+    if (configuration.DoConstrainedMD) {
       // calculate forces and potential
       atom **PermutationMap = NULL;
-      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap, DoConstrained, 0, IsAngstroem);
-      EvaluateConstrainedForces(out, DoConstrained, 0, PermutationMap, &Force);
+      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
+      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
       Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     }
     
     // and perform Verlet integration for each atom with position, velocity and force vector
-    runner = elemente->start;
-    while (runner->next != elemente->end) { // go through every element
-      runner = runner->next;
-      IonMass = runner->mass;
-      a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-      if (ElementsInMolecule[runner->Z]) { // if this element got atoms
-        AtomNo = 0;
+    walker = start;
+    while (walker->next != end) { // go through every atom of this element
+      walker = walker->next;
+      a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
+      // check size of vectors
+      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
+        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
+        Trajectories[walker].R.resize(MDSteps+10);
+        Trajectories[walker].U.resize(MDSteps+10);
+        Trajectories[walker].F.resize(MDSteps+10);
+      }
+      
+      // Update R (and F)
+      for (int d=0; d<NDIM; d++) {
+        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
+        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);
+        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+      }
+      // Update U
+      for (int d=0; d<NDIM; d++) {
+        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
+        Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
+      }
+//      out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
+//      for (int d=0;d<NDIM;d++)
+//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
+//      out << ")" << endl;
+            // next atom
+    }
+  }
+  // correct velocities (rather momenta) so that center of mass remains motionless
+  for(int d=0;d<NDIM;d++)
+    Vector[d] = 0.;
+  IonMass = 0.;
+  walker = start;
+  while (walker->next != end) { // go through every atom
+    walker = walker->next;
+    IonMass += walker->type->mass;  // sum up total mass
+    for(int d=0;d<NDIM;d++) {
+      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
+    }
+  }
+  walker = start;
+  while (walker->next != end) { // go through every atom of this element
+    walker = walker->next;
+    for(int d=0;d<NDIM;d++) {
+      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
+      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
+    }
+  }
+  Thermostats(configuration, ActualTemp, Berendsen);
+  MDSteps++;
+
+  
+  // exit
+  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;
+  double sigma;
+  double IonMass;
+  int d;
+  gsl_rng * r;
+  const gsl_rng_type * T;
+  double *U = NULL, *F = NULL, FConstraint[NDIM];
+  atom *walker = NULL;
+  
+  // calculate scale configuration
+  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
+    
+  // differentating between the various thermostats
+  switch(Thermostat) {
+     case None:
+      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
+      break;
+     case Woodcock:
+      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
+        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
         walker = start;
         while (walker->next != end) { // go through every atom of this element
           walker = walker->next;
-          if (walker->type == runner) { // if this atom fits to element
-            // check size of vectors
-            if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
-              //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-              Trajectories[walker].R.resize(MDSteps+10);
-              Trajectories[walker].U.resize(MDSteps+10);
-              Trajectories[walker].F.resize(MDSteps+10);
+          IonMass = walker->type->mass;
+          U = Trajectories[walker].U.at(MDSteps).x;
+          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              U[d] *= sqrt(ScaleTempFactor);
+              ekin += 0.5*IonMass * U[d]*U[d];
             }
-            
-            // Update R (and F)
-            for (int d=0; d<NDIM; d++) {
-              Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
-              Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
-              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
-              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+        }
+      }
+      break;
+     case Gaussian:
+      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            G += U[d] * F[d]; 
+            E += U[d]*U[d]*IonMass;
+          }
+      }
+      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            FConstraint[d] = (G/E) * (U[d]*IonMass);
+            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
+            ekin += IonMass * U[d]*U[d];
+          }
+      }
+      break;
+     case Langevin:
+      cout << 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
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->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.)) {  
+            cout << Verbose(3) << "Particle " << *walker << " 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 (d=0; d<NDIM; d++) {
+              U[d] = gsl_ran_gaussian (r, sigma);
             }
-            // Update U
-            for (int d=0; d<NDIM; d++) {
-              Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
-              Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
-            }
-//            out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
-//            for (int d=0;d<NDIM;d++)
-//              out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
-//            out << ")\t(";
-//            for (int d=0;d<NDIM;d++)
-//              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
-//            out << ")" << endl;
-            // next atom
-            AtomNo++;
+            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
+          }
+          for (d=0; d<NDIM; d++)
+            ekin += 0.5*IonMass * U[d]*U[d];
+        }
+      }
+      break;
+     case Berendsen:
+      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
+            ekin += 0.5*IonMass * U[d]*U[d];
           }
         }
       }
-    }
-  }
-//  // correct velocities (rather momenta) so that center of mass remains motionless
-//  for(int d=0;d<NDIM;d++)
-//    Vector[d] = 0.;
-//  IonMass = 0.;
-//  walker = start;
-//  while (walker->next != end) { // go through every atom
-//    walker = walker->next;
-//    IonMass += walker->type->mass;  // sum up total mass
-//    for(int d=0;d<NDIM;d++) {
-//      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
-//    }
-//  }
-//  walker = start;
-//  while (walker->next != end) { // go through every atom of this element
-//    walker = walker->next;
-//    for(int d=0;d<NDIM;d++) {
-//      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
-//    }
-//  }
-  MDSteps++;
-
-  
-  // exit
-  return true;
+      break;
+     case NoseHoover:
+      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
+      // dynamically evolve alpha (the additional degree of freedom)
+      delta_alpha = 0.;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            delta_alpha += U[d]*U[d]*IonMass;
+          }
+        }
+      }
+      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
+      configuration.alpha += delta_alpha*configuration.Deltat;
+      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
+      // apply updated alpha as additional force
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
+              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
+              ekin += (0.5*IonMass) * U[d]*U[d];
+            }
+        }
+      }
+      break;
+  }    
+  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
 };
 
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision 16a52b3c50fff2ccc678808e2588e1174646b8ec)
+++ src/molecules.hpp	(revision 62f793aff2f24a09fa44578b18eca5a99c739605)
@@ -16,4 +16,5 @@
 #include <gsl/gsl_multimin.h>
 #include <gsl/gsl_vector.h>
+#include <gsl/gsl_randist.h>
 
 // STL headers
@@ -192,7 +193,13 @@
 };
 
+
 ostream & operator << (ostream &ost, bond &b);
 
 class MoleculeLeafClass;
+
+
+#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.
@@ -260,6 +267,9 @@
 	void PrincipalAxisSystem(ofstream *out, bool DoRotate);
 	double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
-	bool VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained);
-	double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
+	
+	bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
+  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
+
+  double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
 	double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
 	void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
@@ -392,4 +402,12 @@
     int DoConstrainedMD;
     int MaxOuterStep;
+    int Thermostat;
+    int *ThermostatImplemented;
+    char **ThermostatNames;
+    double TempFrequency;
+    double alpha;
+    double HooverMass;
+    double TargetTemp;
+    int ScaleTempStep;
     
     private:
@@ -415,6 +433,4 @@
     int OutVisStep;
     int OutSrcStep;
-    double TargetTemp;
-    int ScaleTempStep;
     int MaxPsiStep;
     double EpsWannier;
@@ -463,4 +479,5 @@
   char *GetDefaultPath() const;
   void SetDefaultPath(const char *path);
+  void InitThermostats(ifstream *source);
 };
 
