Index: src/Thermostats/Berendsen.cpp
===================================================================
--- src/Thermostats/Berendsen.cpp	(revision a5028f9fe5de8314faeb80f8192d465807b8d3f0)
+++ src/Thermostats/Berendsen.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -63,6 +63,6 @@
   double ScaleTempFactor = getContainer().TargetTemp/ActualTemp;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       U *= sqrt(1+(World::getInstance().getConfig()->Deltat/TempFrequency)*(ScaleTempFactor-1));
       ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
Index: src/Thermostats/GaussianThermostat.cpp
===================================================================
--- src/Thermostats/GaussianThermostat.cpp	(revision a5028f9fe5de8314faeb80f8192d465807b8d3f0)
+++ src/Thermostats/GaussianThermostat.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -70,6 +70,6 @@
   double ekin =0;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) {// even FixedIon moves, only not by other's forces
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) {// even FixedIon moves, only not by other's forces
       U += World::getInstance().getConfig()->Deltat * G_over_E * U;
       ekin += (*iter)->getType()->getMass() * U.NormSquared();
@@ -84,7 +84,7 @@
   G=0;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    Vector &F = (*iter)->Trajectory.F.at(step);
-    if ((*iter)->FixedIon == 0){ // even FixedIon moves, only not by other's forces
+    const Vector &U = (*iter)->getAtomicVelocity(step);
+    const Vector &F = (*iter)->getAtomicForce(step);
+    if ((*iter)->getFixedIon() == 0){ // even FixedIon moves, only not by other's forces
       G += U.ScalarProduct(F);
       E += U.NormSquared()*(*iter)->getType()->getMass();
Index: src/Thermostats/Langevin.cpp
===================================================================
--- src/Thermostats/Langevin.cpp	(revision a5028f9fe5de8314faeb80f8192d465807b8d3f0)
+++ src/Thermostats/Langevin.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -20,4 +20,8 @@
 #include "ThermoStatContainer.hpp"
 
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+
 #include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
 #include "RandomNumbers/RandomNumberGenerator.hpp"
@@ -25,19 +29,20 @@
 Langevin::Langevin(double _TempFrequency,double _alpha) :
   TempFrequency(_TempFrequency),
-  alpha(_alpha)
-{
-  gsl_rng_env_setup();
-  T = gsl_rng_default;
-  r = gsl_rng_alloc (T);
-}
+  alpha(_alpha),
+  rng_engine(new boost::mt19937),
+  rng_distribution(NULL)
+{}
 
 Langevin::Langevin() :
   TempFrequency(2.5),
-  alpha(0.)
+  alpha(0.),
+  rng_engine(new boost::mt19937),
+  rng_distribution(NULL)
 {}
 
 Langevin::~Langevin()
 {
-  gsl_rng_free (r);
+  delete rng_engine;
+  delete rng_distribution;
 }
 
@@ -76,5 +81,5 @@
 double Langevin::doScaleAtoms(unsigned int step,double ActualTemp,ForwardIterator begin, ForwardIterator end){
   DoLog(2) && (Log() << Verbose(2) <<  "Applying Langevin thermostat..." << endl);
-  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
+  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator("mt19937", "normal_distribution");
   const double rng_min = random.min();
   const double rng_max = random.max();
@@ -82,6 +87,8 @@
   for(ForwardIterator iter=begin;iter!=end;++iter){
     double sigma  = sqrt(getContainer().TargetTemp/(*iter)->getType()->getMass()); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    rng_distribution = new boost::normal_distribution<>(0,sigma);
+    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > rng(*rng_engine, *rng_distribution);
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 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 (((((random()/(rng_max-rng_min)))*TempFrequency) < 1.)) {
@@ -89,5 +96,5 @@
         // 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);
+          U[d] = rng();
         }
         DoLog(2) && (Log() << Verbose(2) << U.Norm() << endl);
@@ -95,4 +102,6 @@
       ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
     }
+    delete rng_distribution;
+    rng_distribution = NULL;
   }
   return ekin;
Index: src/Thermostats/Langevin.hpp
===================================================================
--- src/Thermostats/Langevin.hpp	(revision a5028f9fe5de8314faeb80f8192d465807b8d3f0)
+++ src/Thermostats/Langevin.hpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -14,6 +14,8 @@
 #endif
 
+#include "Thermostats/Thermostat.hpp"
 
-#include "Thermostats/Thermostat.hpp"
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
 
 class Langevin : public Thermostat
@@ -38,6 +40,7 @@
   double TempFrequency;
   double alpha;
-  gsl_rng * r;
-  const gsl_rng_type * T;
+
+  boost::mt19937 *rng_engine;
+  boost::normal_distribution<> *rng_distribution;
 };
 
Index: src/Thermostats/NoseHoover.cpp
===================================================================
--- src/Thermostats/NoseHoover.cpp	(revision a5028f9fe5de8314faeb80f8192d465807b8d3f0)
+++ src/Thermostats/NoseHoover.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -68,6 +68,6 @@
   double ekin =0;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       U += World::getInstance().getConfig()->Deltat/(*iter)->getType()->getMass() * (alpha * (U * (*iter)->getType()->getMass()));
       ekin += (0.5*(*iter)->getType()->getMass()) * U.NormSquared();
@@ -82,6 +82,6 @@
   count=0;
   for(ForwardIterator iter = begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    const Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       delta_alpha += U.NormSquared()*(*iter)->getType()->getMass();
     }
Index: src/Thermostats/Woodcock.cpp
===================================================================
--- src/Thermostats/Woodcock.cpp	(revision a5028f9fe5de8314faeb80f8192d465807b8d3f0)
+++ src/Thermostats/Woodcock.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -67,6 +67,6 @@
     double ekin;
     for (ForwardIterator iter = begin; iter!=end;++iter){
-      Vector &U = (*iter)->Trajectory.U.at(step);
-      if ((*iter)->FixedIon == 0){ // even FixedIon moves, only not by other's forces
+      Vector &U = (*iter)->getAtomicVelocity(step);
+      if ((*iter)->getFixedIon() == 0){ // even FixedIon moves, only not by other's forces
         U *= ScaleTempFactor;
         ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
