Index: src/AtomSet.hpp
===================================================================
--- src/AtomSet.hpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/AtomSet.hpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -82,9 +82,9 @@
   template<class T>
   struct valueSum {
-    valueSum(T (atom::*_f)() const,T startValue) :
+    valueSum(T (AtomInfo::*_f)() const,T startValue) :
       f(_f),
       value(startValue)
     {}
-    T operator+(atom *atom){
+    T operator+(AtomInfo *atom){
       return value + (atom->*f)();
     }
@@ -93,5 +93,5 @@
       return value;
     }
-    T (atom::*f)() const;
+    T (AtomInfo::*f)() const;
     T value;
   };
@@ -99,10 +99,10 @@
   template<class T>
   struct stepValueSum {
-    stepValueSum(unsigned int _step, T (atom::*_f)(unsigned int) const,T startValue) :
+    stepValueSum(unsigned int _step, T (AtomInfo::*_f)(unsigned int) const,T startValue) :
       step(_step),
       f(_f),
       value(startValue)
     {}
-    T operator+(atom *atom){
+    T operator+(AtomInfo *atom){
       return value + (atom->*f)(step);
     }
@@ -112,5 +112,5 @@
     }
     unsigned int step;
-    T (atom::*f)(unsigned int) const;
+    T (AtomInfo::*f)(unsigned int) const;
     T value;
   };
@@ -119,5 +119,5 @@
 template<class Set>
 inline void AtomSetMixin<Set>::translate(const Vector &translater){
-  BOOST_FOREACH(atom *atom,*this){
+  BOOST_FOREACH(AtomInfo *atom,*this){
     *(atom) += translater;
   }
@@ -126,6 +126,6 @@
 template<class Set>
 inline void AtomSetMixin<Set>::addVelocityAtStep(const Vector velocity, unsigned int step){
-  BOOST_FOREACH(atom *atom,*this){
-    atom->Trajectory.U.at(step) += velocity;
+  BOOST_FOREACH(AtomInfo *atom,*this){
+    atom->getAtomicVelocity(step) += velocity;
   }
 }
@@ -141,15 +141,15 @@
 template<class Set>
 inline double AtomSetMixin<Set>::totalMass() const{
-  return accumulate(this->begin(),this->end(),valueSum<double>(&atom::getMass,0)).value;
+  return accumulate(this->begin(),this->end(),valueSum<double>(&AtomInfo::getMass,0)).value;
 }
 
 template<class Set>
 inline double AtomSetMixin<Set>::totalTemperatureAtStep(unsigned int step) const{
-  return accumulate(this->begin(),this->end(),stepValueSum<double>(step,&atom::getKineticEnergy,0)).value;
+  return accumulate(this->begin(),this->end(),stepValueSum<double>(step,&AtomInfo::getKineticEnergy,0)).value;
 }
 
 template<class Set>
 inline Vector AtomSetMixin<Set>::totalMomentumAtStep(unsigned int step) const{
-  return accumulate(this->begin(),this->end(),stepValueSum<Vector>(step,&atom::getMomentum,Vector())).value;
+  return accumulate(this->begin(),this->end(),stepValueSum<Vector>(step,&AtomInfo::getMomentum,Vector())).value;
 }
 
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/Makefile.am	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -21,7 +21,5 @@
   atom_graphnode.cpp \
   atom_graphnodeinfo.cpp \
-  atom_particleinfo.cpp \
-  atom_trajectoryparticle.cpp \
-  atom_trajectoryparticleinfo.cpp
+  atom_particleinfo.cpp
 ATOMHEADER = \
   atom.hpp \
@@ -32,7 +30,5 @@
   atom_graphnode.hpp \
   atom_graphnodeinfo.hpp \
-  atom_particleinfo.hpp \
-  atom_trajectoryparticle.hpp \
-  atom_trajectoryparticleinfo.hpp
+  atom_particleinfo.hpp
 
 ANALYSISSOURCE = \
Index: src/Parser/PcpParser.cpp
===================================================================
--- src/Parser/PcpParser.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/Parser/PcpParser.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -557,5 +557,5 @@
     *file << "Ion_Type" << ZtoIndexMap[Z] << "_" << ZtoCountMap[Z] << "\t"  << fixed << setprecision(9) << showpoint;
     *file << (*AtomRunner)->at(0) << "\t" << (*AtomRunner)->at(1) << "\t" << (*AtomRunner)->at(2);
-    *file << "\t" << (*AtomRunner)->FixedIon;
+    *file << "\t" << (int)((*AtomRunner)->getFixedIon());
     if ((*AtomRunner)->getAtomicVelocity().Norm() > MYEPSILON)
       *file << "\t" << scientific << setprecision(6) << (*AtomRunner)->getAtomicVelocity()[0] << "\t" << (*AtomRunner)->getAtomicVelocity()[1] << "\t" << (*AtomRunner)->getAtomicVelocity()[2] << "\t";
Index: src/Thermostats/Berendsen.cpp
===================================================================
--- src/Thermostats/Berendsen.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 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 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 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 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 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 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 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 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 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 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 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();
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/atom.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -55,5 +55,5 @@
   setPosition(pointer->getPosition()); // copy coordination
   setAtomicVelocity(pointer->getAtomicVelocity()); // copy velocity
-  FixedIon = pointer->FixedIon;
+  setFixedIon(pointer->getFixedIon());
 };
 
@@ -152,5 +152,5 @@
     *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
     *out << at(0) << "\t" << at(1) << "\t" << at(2);
-    *out << "\t" << FixedIon;
+    *out << "\t" << (int)(getFixedIon());
     if (getAtomicVelocity().Norm() > MYEPSILON)
       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
@@ -179,5 +179,5 @@
     *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
     *out << at(0) << "\t" << at(1) << "\t" << at(2);
-    *out << "\t" << FixedIon;
+    *out << "\t" << getFixedIon();
     if (getAtomicVelocity().Norm() > MYEPSILON)
       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
@@ -191,5 +191,5 @@
 };
 
-/** Output of a single atom as one lin in xyz file.
+/** Output of a single atom as one line in xyz file.
  * \param *out stream to output to
   * \return true - \a *out present, false - \a *out is NULL
@@ -204,5 +204,5 @@
 };
 
-/** Output of a single atom as one lin in xyz file.
+/** Output of a single atom as one line in xyz file.
  * \param *out stream to output to
  * \param *ElementNo array with ion type number in the config file this atom's element shall have
@@ -218,10 +218,10 @@
     ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration");
     *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[getType()->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
-    *out << Trajectory.R.at(step)[0] << "\t" << Trajectory.R.at(step)[1] << "\t" << Trajectory.R.at(step)[2];
-    *out << "\t" << FixedIon;
-    if (Trajectory.U.at(step).Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step)[0] << "\t" << Trajectory.U.at(step)[1] << "\t" << Trajectory.U.at(step)[2] << "\t";
-    if (Trajectory.F.at(step).Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step)[0] << "\t" << Trajectory.F.at(step)[1] << "\t" << Trajectory.F.at(step)[2] << "\t";
+    *out << getPosition(step)[0] << "\t" << getPosition(step)[1] << "\t" << getPosition(step)[2];
+    *out << "\t" << (int)(getFixedIon());
+    if (getAtomicVelocity(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << getAtomicVelocity(step)[0] << "\t" << getAtomicVelocity(step)[1] << "\t" << getAtomicVelocity(step)[2] << "\t";
+    if (getAtomicForce(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << getAtomicForce(step)[0] << "\t" << getAtomicForce(step)[1] << "\t" << getAtomicForce(step)[2] << "\t";
     *out << "\t# Number in molecule " << nr << endl;
     return true;
@@ -239,7 +239,7 @@
   if (out != NULL) {
     *out << getType()->getSymbol() << "\t";
-    *out << Trajectory.R.at(step)[0] << "\t";
-    *out << Trajectory.R.at(step)[1] << "\t";
-    *out << Trajectory.R.at(step)[2] << endl;
+    *out << getPosition(step)[0] << "\t";
+    *out << getPosition(step)[1] << "\t";
+    *out << getPosition(step)[2] << endl;
     return true;
   } else
@@ -383,8 +383,4 @@
 }
 
-double atom::getMass() const{
-  return getType()->getMass();
-}
-
 int atom::getNr() const{
   return nr;
Index: src/atom.hpp
===================================================================
--- src/atom.hpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/atom.hpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -27,5 +27,4 @@
 #include "atom_graphnode.hpp"
 #include "atom_particleinfo.hpp"
-#include "atom_trajectoryparticle.hpp"
 #include "TesselPoint.hpp"
 #include "types.hpp"
@@ -45,5 +44,5 @@
  * Class incorporates position, type
  */
-class atom : public TrajectoryParticle, public GraphNode, public BondedParticle, public TesselPoint {
+class atom : public GraphNode, public BondedParticle, public TesselPoint {
   friend atom* NewAtom(atomId_t);
   friend void  DeleteAtom(atom*);
@@ -98,5 +97,4 @@
    void removeFromMolecule();
 
-   double getMass() const;
    int getNr() const;
 
Index: src/atom_atominfo.cpp
===================================================================
--- src/atom_atominfo.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/atom_atominfo.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -20,7 +20,11 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "config.hpp"
+#include "element.hpp"
+#include "parser.hpp"
 #include "periodentafel.hpp"
 #include "World.hpp"
-#include "element.hpp"
 #include "atom_atominfo.hpp"
 
@@ -28,5 +32,6 @@
  */
 AtomInfo::AtomInfo() :
-  AtomicElement(NULL)
+  AtomicElement(NULL),
+  FixedIon(false)
 {
   AtomicPosition.reserve(1);
@@ -42,11 +47,22 @@
 AtomInfo::AtomInfo(const AtomInfo &_atom) :
     AtomicPosition(_atom.AtomicPosition),
-    AtomicElement(_atom.AtomicElement)
-{};
+    AtomicElement(_atom.AtomicElement),
+    FixedIon(false)
+{
+  AtomicVelocity.reserve(1);
+  AtomicVelocity.push_back(zeroVec);
+  AtomicForce.reserve(1);
+  AtomicForce.push_back(zeroVec);
+};
 
 AtomInfo::AtomInfo(const VectorInterface &_v) :
-    AtomicElement(NULL)
+    AtomicElement(NULL),
+    FixedIon(false)
 {
   AtomicPosition[0] = _v.getPosition();
+  AtomicVelocity.reserve(1);
+  AtomicVelocity.push_back(zeroVec);
+  AtomicForce.reserve(1);
+  AtomicForce.push_back(zeroVec);
 };
 
@@ -80,4 +96,11 @@
 {
   return AtomicPosition[0];
+}
+
+const Vector& AtomInfo::getPosition(const int _step) const
+{
+  ASSERT(_step < AtomicPosition.size(),
+      "AtomInfo::getPosition() - Access out of range.");
+  return AtomicPosition[_step];
 }
 
@@ -96,4 +119,11 @@
 }
 
+Vector& AtomInfo::getAtomicVelocity(const int _step)
+{
+  ASSERT(_step < AtomicVelocity.size(),
+      "AtomInfo::getAtomicVelocity() - Access out of range.");
+  return AtomicVelocity[_step];
+}
+
 const Vector& AtomInfo::getAtomicVelocity() const
 {
@@ -101,4 +131,11 @@
 }
 
+const Vector& AtomInfo::getAtomicVelocity(const int _step) const
+{
+  ASSERT(_step < AtomicVelocity.size(),
+      "AtomInfo::getAtomicVelocity() - Access out of range.");
+  return AtomicVelocity[_step];
+}
+
 void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
 {
@@ -106,4 +143,15 @@
 }
 
+void AtomInfo::setAtomicVelocity(const int _step, const Vector &_newvelocity)
+{
+  ASSERT(_step <= AtomicVelocity.size(),
+      "AtomInfo::setAtomicVelocity() - Access out of range.");
+  if(_step < (int)AtomicVelocity.size()) {
+    AtomicVelocity[_step] = _newvelocity;
+  } else if (_step == (int)AtomicVelocity.size()) {
+    AtomicVelocity.push_back(_newvelocity);
+  }
+}
+
 const Vector& AtomInfo::getAtomicForce() const
 {
@@ -111,4 +159,11 @@
 }
 
+const Vector& AtomInfo::getAtomicForce(const int _step) const
+{
+  ASSERT(_step < AtomicForce.size(),
+      "AtomInfo::getAtomicForce() - Access out of range.");
+  return AtomicForce[_step];
+}
+
 void AtomInfo::setAtomicForce(const Vector &_newforce)
 {
@@ -116,7 +171,40 @@
 }
 
+void AtomInfo::setAtomicForce(const int _step, const Vector &_newforce)
+{
+  ASSERT(_step <= AtomicForce.size(),
+      "AtomInfo::setAtomicForce() - Access out of range.");
+  if(_step < (int)AtomicForce.size()) {
+    AtomicForce[_step] = _newforce;
+  } else if (_step == (int)AtomicForce.size()) {
+    AtomicForce.push_back(_newforce);
+  }
+}
+
+bool AtomInfo::getFixedIon() const
+{
+  return FixedIon;
+}
+
+void AtomInfo::setFixedIon(const bool _fixedion)
+{
+  FixedIon = _fixedion;
+}
+
 void AtomInfo::setPosition(const Vector& _vector)
 {
   AtomicPosition[0] = _vector;
+  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
+}
+
+void AtomInfo::setPosition(int _step, const Vector& _vector)
+{
+  ASSERT(_step <= AtomicPosition.size(),
+      "AtomInfo::setPosition() - Access out of range.");
+  if(_step < (int)AtomicPosition.size()) {
+    AtomicPosition[_step] = _vector;
+  } else if (_step == (int)AtomicPosition.size()) {
+    AtomicPosition.push_back(_vector);
+  }
   //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
 }
@@ -204,4 +292,90 @@
 }
 
+/**
+ *  returns the kinetic energy of this atom at a given time step
+ */
+double AtomInfo::getKineticEnergy(unsigned int step) const{
+  return getMass() * AtomicVelocity[step].NormSquared();
+}
+
+Vector AtomInfo::getMomentum(unsigned int step) const{
+  return getMass()*AtomicVelocity[step];
+}
+
+/** Extends the trajectory STL vector to the new size.
+ * Does nothing if \a MaxSteps is smaller than current size.
+ * \param MaxSteps
+ */
+void AtomInfo::ResizeTrajectory(size_t MaxSteps)
+{
+  if (AtomicPosition.size() <= (unsigned int)(MaxSteps)) {
+    DoLog(0) && (Log() << Verbose(0) << "Increasing size for trajectory array from " << AtomicPosition.size() << " to " << (MaxSteps+1) << "." << endl);
+    AtomicPosition.resize(MaxSteps+1, zeroVec);
+    AtomicVelocity.resize(MaxSteps+1, zeroVec);
+    AtomicForce.resize(MaxSteps+1, zeroVec);
+  }
+};
+
+size_t AtomInfo::getTrajectorySize() const
+{
+  return AtomicPosition.size();
+}
+
+double AtomInfo::getMass() const{
+  return AtomicElement->getMass();
+}
+
+/** Copies a given trajectory step \a src onto another \a dest
+ * \param dest index of destination step
+ * \param src index of source step
+ */
+void AtomInfo::CopyStepOnStep(int dest, int src)
+{
+  if (dest == src)  // self assignment check
+    return;
+
+  ASSERT(dest < AtomicPosition.size(),
+      "AtomInfo::CopyStepOnStep() - destination outside of current trajectory array.");
+  ASSERT(src < AtomicPosition.size(),
+      "AtomInfo::CopyStepOnStep() - source outside of current trajectory array.");
+  for (int n=NDIM;n--;) {
+    AtomicPosition.at(dest)[n] = AtomicPosition.at(src)[n];
+    AtomicVelocity.at(dest)[n] = AtomicVelocity.at(src)[n];
+    AtomicForce.at(dest)[n] = AtomicForce.at(src)[n];
+  }
+};
+
+/** Performs a velocity verlet update of the trajectory.
+ * Parameters are according to those in configuration class.
+ * \param NextStep index of sequential step to set
+ * \param *configuration pointer to configuration with parameters
+ * \param *Force matrix with forces
+ */
+void AtomInfo::VelocityVerletUpdate(int nr, int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
+{
+  //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
+  if ((int)AtomicPosition.size() <= NextStep)
+    ResizeTrajectory(NextStep+10);
+  for (int d=0; d<NDIM; d++) {
+    AtomicForce.at(NextStep)[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    AtomicPosition.at(NextStep)[d] = AtomicPosition.at(NextStep-1)[d];
+    AtomicPosition.at(NextStep)[d] += configuration->Deltat*(AtomicVelocity.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+    AtomicPosition.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(AtomicForce.at(NextStep)[d]/getMass());     // F = m * a and s =
+  }
+  // Update U
+  for (int d=0; d<NDIM; d++) {
+    AtomicVelocity.at(NextStep)[d] = AtomicVelocity.at(NextStep-1)[d];
+    AtomicVelocity.at(NextStep)[d] += configuration->Deltat * (AtomicForce.at(NextStep)[d]+AtomicForce.at(NextStep-1)[d]/getMass()); // v = F/m * t
+  }
+  // Update R (and F)
+//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
+//      for (int d=0;d<NDIM;d++)
+//        out << AtomicPosition.at(NextStep).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        Log() << Verbose(0) << AtomicVelocity.at(NextStep).x[d] << " ";          // next step
+//      out << ")" << endl;
+};
+
 const AtomInfo& operator*=(AtomInfo& a, const double m)
 {
Index: src/atom_atominfo.hpp
===================================================================
--- src/atom_atominfo.hpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/atom_atominfo.hpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -27,5 +27,7 @@
 
 class AtomInfo;
+class config;
 class element;
+class ForceMatrix;
 class RealSpaceMatrix;
 
@@ -58,4 +60,6 @@
   /** Getter for AtomicVelocity.
    *
+   * Current time step is used.
+   *
    * @return constant reference to AtomicVelocity
    */
@@ -63,34 +67,146 @@
   /** Getter for AtomicVelocity.
    *
+   * @param _step time step to return
+   * @return constant reference to AtomicVelocity
+   */
+  Vector& getAtomicVelocity(const int _step);
+  /** Getter for AtomicVelocity.
+   *
+   * Current time step is used.
+   *
    * @return constant reference to AtomicVelocity
    */
   const Vector& getAtomicVelocity() const;
+  /** Getter for AtomicVelocity.
+   *
+   * @param _step time step to return
+   * @return constant reference to AtomicVelocity
+   */
+  const Vector& getAtomicVelocity(const int _step) const;
   /** Setter for AtomicVelocity.
    *
+   * Current time step is used.
+   *
    * @param _newvelocity new velocity to set
    */
   void setAtomicVelocity(const Vector &_newvelocity);
+  /** Setter for AtomicVelocity.
+   *
+   * @param _step time step to set
+   * @param _newvelocity new velocity to set
+   */
+  void setAtomicVelocity(const int _step, const Vector &_newvelocity);
 
   /** Getter for AtomicForce.
    *
+   * Current time step is used.
+   *
    * @return constant reference to AtomicForce
    */
   const Vector& getAtomicForce() const;
+  /** Getter for AtomicForce.
+   *
+   * @param _step time step to return
+   * @return constant reference to AtomicForce
+   */
+  const Vector& getAtomicForce(const int _step) const;
   /** Setter for AtomicForce.
    *
+   * Current time step is used.
+   *
    * @param _newvelocity new force vector to set
    */
   void setAtomicForce(const Vector &_newforce);
+  /** Setter for AtomicForce.
+   *
+   * @param _step time step to set
+   * @param _newvelocity new force vector to set
+   */
+  void setAtomicForce(const int _step, const Vector &_newforce);
+
+  /** Getter for FixedIon.
+   *
+   * @return constant reference to FixedIon
+   */
+  bool getFixedIon() const;
+  /** Setter for FixedIon.
+   *
+   * @param _fixedion new state of FixedIon
+   */
+  void setFixedIon(const bool _fixedion);
 
   ///// manipulation of the atomic position
 
   // Accessors ussually come in pairs... and sometimes even more than that
+  /** Getter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @param i component of vector
+   * @return i-th component of atomic position
+   */
   const double& operator[](size_t i) const;
+  /** Getter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * \sa operator[], this is if instance is a reference.
+   *
+   * @param i component of vector
+   * @return i-th component of atomic position
+   */
   const double& at(size_t i) const;
+  /** Getter for AtomicPosition.
+   *
+   * \sa operator[], this is if instance is a reference.
+   *
+   * @param i index of component of AtomicPosition
+   * @param _step time step to return
+   * @return atomic position at time step _step
+   */
+  const Vector& atStep(size_t i, int _step) const;
+  /** Setter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @param i component to set
+   * @param value value to set to
+   */
   void set(size_t i, const double value);
+  /** Setter for AtomicPosition.
+   *
+   * @param i component to set
+   * @param _step time step to set
+   * @param value value to set to
+   */
+  void setAtStep(size_t i, int _step, const double value);
+  /** Getter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @return atomic position
+   */
   const Vector& getPosition() const;
+  /** Getter for AtomicPosition.
+   *
+   * @param _step time step to return
+   * @return atomic position at time step _step
+   */
+  const Vector& getPosition(int _step) const;
 
   // Assignment operator
+  /** Setter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @param _vector new position to set
+   */
   void setPosition(const Vector& _vector);
+  /** Setter for AtomicPosition.
+   *
+   * @param _step time step to set
+   * @param _vector new position to set for time step _step
+   */
+  void setPosition(int _step, const Vector& _vector);
   class VectorInterface &operator=(const Vector& _vector);
 
@@ -114,4 +230,13 @@
   void Scale(const double factor);
 
+  // operations for trajectories
+  void ResizeTrajectory(size_t MaxSteps);
+  size_t getTrajectorySize() const;
+  void CopyStepOnStep(int dest, int src);
+  void VelocityVerletUpdate(int nr, int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
+  double getKineticEnergy(unsigned int step) const;
+  Vector getMomentum(unsigned int step) const;
+  double getMass() const;
+
   std::ostream & operator << (std::ostream &ost) const;
 
@@ -122,4 +247,5 @@
 
   const element *AtomicElement;  //!< pointing to element
+  bool FixedIon;
 };
 
Index: src/atom_trajectoryparticle.cpp
===================================================================
--- src/atom_trajectoryparticle.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 	(revision )
@@ -1,129 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010 University of Bonn. All rights reserved.
- * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
- */
-
-/*
- * atom_trajectoryparticle.cpp
- *
- *  Created on: Oct 19, 2009
- *      Author: heber
- */
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "CodePatterns/MemDebug.hpp"
-
-#include "atom.hpp"
-#include "atom_trajectoryparticle.hpp"
-#include "config.hpp"
-#include "element.hpp"
-#include "CodePatterns/Info.hpp"
-#include "CodePatterns/Log.hpp"
-#include "parser.hpp"
-#include "ThermoStatContainer.hpp"
-#include "CodePatterns/Verbose.hpp"
-
-/** Constructor of class TrajectoryParticle.
- */
-TrajectoryParticle::TrajectoryParticle()
-{
-};
-
-/** Destructor of class TrajectoryParticle.
- */
-TrajectoryParticle::~TrajectoryParticle()
-{
-};
-
-/**
- *  returns the kinetic energy of this atom at a given time step
- */
-double TrajectoryParticle::getKineticEnergy(unsigned int step) const{
-  return getType()->getMass() * Trajectory.U.at(step).NormSquared();
-}
-
-Vector TrajectoryParticle::getMomentum(unsigned int step) const{
-  return getType()->getMass()*Trajectory.U.at(step);
-}
-
-/** Extends the trajectory STL vector to the new size.
- * Does nothing if \a MaxSteps is smaller than current size.
- * \param MaxSteps
- */
-void TrajectoryParticle::ResizeTrajectory(int MaxSteps)
-{
-  Info FunctionInfo(__func__);
-  if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
-    DoLog(0) && (Log() << Verbose(0) << "Increasing size for trajectory array of " << nr << " from " << Trajectory.R.size() << " to " << (MaxSteps+1) << "." << endl);
-    Trajectory.R.resize(MaxSteps+1);
-    Trajectory.U.resize(MaxSteps+1);
-    Trajectory.F.resize(MaxSteps+1);
-  }
-};
-
-/** Copies a given trajectory step \a src onto another \a dest
- * \param dest index of destination step
- * \param src index of source step
- */
-void TrajectoryParticle::CopyStepOnStep(int dest, int src)
-{
-  if (dest == src)  // self assignment check
-    return;
-
-  for (int n=NDIM;n--;) {
-    Trajectory.R.at(dest)[n] = Trajectory.R.at(src)[n];
-    Trajectory.U.at(dest)[n] = Trajectory.U.at(src)[n];
-    Trajectory.F.at(dest)[n] = Trajectory.F.at(src)[n];
-  }
-};
-
-/** Performs a velocity verlet update of the trajectory.
- * Parameters are according to those in configuration class.
- * \param NextStep index of sequential step to set
- * \param *configuration pointer to configuration with parameters
- * \param *Force matrix with forces
- */
-void TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
-{
-  //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
-  for (int d=0; d<NDIM; d++) {
-    Trajectory.F.at(NextStep)[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-    Trajectory.R.at(NextStep)[d] = Trajectory.R.at(NextStep-1)[d];
-    Trajectory.R.at(NextStep)[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-    Trajectory.R.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep)[d]/getType()->getMass());     // F = m * a and s =
-  }
-  // Update U
-  for (int d=0; d<NDIM; d++) {
-    Trajectory.U.at(NextStep)[d] = Trajectory.U.at(NextStep-1)[d];
-    Trajectory.U.at(NextStep)[d] += configuration->Deltat * (Trajectory.F.at(NextStep)[d]+Trajectory.F.at(NextStep-1)[d]/getType()->getMass()); // v = F/m * t
-  }
-  // Update R (and F)
-//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
-//      for (int d=0;d<NDIM;d++)
-//        out << Trajectory.R.at(NextStep).x[d] << " ";          // next step
-//      out << ")\t(";
-//      for (int d=0;d<NDIM;d++)
-//        Log() << Verbose(0) << Trajectory.U.at(NextStep).x[d] << " ";          // next step
-//      out << ")" << endl;
-};
-
-std::ostream & TrajectoryParticle::operator << (std::ostream &ost) const
-{
-  ParticleInfo::operator<<(ost);
-  ost << "," << getPosition();
-  return ost;
-}
-
-std::ostream & operator << (std::ostream &ost, const TrajectoryParticle &a)
-{
-  a.ParticleInfo::operator<<(ost);
-  ost << "," << a.getPosition();
-  return ost;
-}
-
Index: src/atom_trajectoryparticle.hpp
===================================================================
--- src/atom_trajectoryparticle.hpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 	(revision )
@@ -1,57 +1,0 @@
-/*
- * atom_trajectoryparticle.hpp
- *
- *  Created on: Oct 19, 2009
- *      Author: heber
- */
-
-#ifndef ATOM_TRAJECTORYPARTICLE_HPP_
-#define ATOM_TRAJECTORYPARTICLE_HPP_
-
-using namespace std;
-
-/*********************************************** includes ***********************************/
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <fstream>
-
-#include <gsl/gsl_randist.h>
-
-#include "atom_atominfo.hpp"
-#include "atom_particleinfo.hpp"
-#include "atom_trajectoryparticleinfo.hpp"
-
-/****************************************** forward declarations *****************************/
-
-class atom;
-class config;
-class ForceMatrix;
-class Vector;
-
-/********************************************** declarations *******************************/
-
-
-class TrajectoryParticle :  public TrajectoryParticleInfo, public virtual AtomInfo, public virtual ParticleInfo {
-public:
-
-  TrajectoryParticle();
-  virtual ~TrajectoryParticle();
-
-  void ResizeTrajectory(int MaxSteps);
-  void CopyStepOnStep(int dest, int src);
-  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
-  double getKineticEnergy(unsigned int step) const;
-  Vector getMomentum(unsigned int step) const;
-
-  std::ostream & operator << (std::ostream &ost) const;
-
-private:
-};
-
-ostream & operator << (ostream &ost, const TrajectoryParticle &a);
-
-#endif /* ATOM_TRAJECTORYPARTICLE_HPP_ */
Index: src/atom_trajectoryparticleinfo.cpp
===================================================================
--- src/atom_trajectoryparticleinfo.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 	(revision )
@@ -1,37 +1,0 @@
-/*
- * Project: MoleCuilder
- * Description: creates and alters molecular systems
- * Copyright (C)  2010 University of Bonn. All rights reserved.
- * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
- */
-
-/*
- * atom_trajectoryparticleinfo.cpp
- *
- *  Created on: Oct 19, 2009
- *      Author: heber
- */
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "CodePatterns/MemDebug.hpp"
-
-#include "atom_trajectoryparticleinfo.hpp"
-
-/** Constructor of class TrajectoryParticleInfo.
- */
-TrajectoryParticleInfo::TrajectoryParticleInfo() :
-  FixedIon(0)
-{};
-
-/** Destructor of class TrajectoryParticleInfo.
- */
-TrajectoryParticleInfo::~TrajectoryParticleInfo()
-{
-  Trajectory.R.clear();
-  Trajectory.U.clear();
-  Trajectory.F.clear();
-};
Index: src/atom_trajectoryparticleinfo.hpp
===================================================================
--- src/atom_trajectoryparticleinfo.hpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ 	(revision )
@@ -1,43 +1,0 @@
-/*
- * atom_trajectoryparticleinfo.hpp
- *
- *  Created on: Oct 19, 2009
- *      Author: heber
- */
-
-#ifndef ATOM_TRAJECTORYPARTICLEINFO_HPP_
-#define ATOM_TRAJECTORYPARTICLEINFO_HPP_
-
-using namespace std;
-
-/*********************************************** includes ***********************************/
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <vector>
-
-#include "LinearAlgebra/Vector.hpp"
-#include "LinearAlgebra/VectorSet.hpp"
-
-/********************************************** declarations *******************************/
-
-class TrajectoryParticleInfo {
-public:
-  struct
-  {
-    VectorSet<vector<Vector> > R;  //!< position vector
-    VectorSet<vector<Vector> > U;  //!< velocity vector
-    VectorSet<vector<Vector> > F;  //!< last force vector
-  } Trajectory;
-  int FixedIon;   //!< config variable that states whether forces act on the ion or not
-
-  TrajectoryParticleInfo();
-  virtual ~TrajectoryParticleInfo();
-
-private:
-};
-
-#endif /* ATOM_TRAJECTORYPARTICLEINFO_HPP_ */
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/config.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -462,4 +462,6 @@
   int No[MAX_ELEMENTS];
   int verbose = 0;
+  Vector tempVector;
+  int _fixedion;
   double value[3];
 
@@ -525,50 +527,42 @@
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], 1, (repetition == 0) ? critical : optional) &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], 1, (repetition == 0) ? critical : optional) &&
-                    ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &_fixedion, 1, (repetition == 0) ? critical : optional));
             if (!status)
               break;
-            neues ->setPosition(position);
+            neues->setFixedIon(_fixedion);
+            neues->setPosition(position);
 
             // check size of vectors
-            if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
+            if (neues->getTrajectorySize() <= (unsigned int)(repetition)) {
               //Log() << Verbose(0) << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
-              neues->Trajectory.R.resize(repetition+10);
-              neues->Trajectory.U.resize(repetition+10);
-              neues->Trajectory.F.resize(repetition+10);
+              neues->ResizeTrajectory(repetition+1);
+              neues->ResizeTrajectory(repetition+1);
+              neues->ResizeTrajectory(repetition+1);
             }
 
-            // put into trajectories list
-            for (int d=0;d<NDIM;d++)
-              neues->Trajectory.R.at(repetition)[d] = neues->at(d);
-
             // parse velocities if present
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->getAtomicVelocity()[0], 1,optional))
-              neues->getAtomicVelocity()[0] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->getAtomicVelocity()[1], 1,optional))
-              neues->getAtomicVelocity()[1] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->getAtomicVelocity()[2], 1,optional))
-              neues->getAtomicVelocity()[2] = 0.;
-            for (int d=0;d<NDIM;d++)
-              neues->Trajectory.U.at(repetition)[d] = neues->getAtomicVelocity()[d];
+            tempVector.Zero();
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &tempVector[0], 1,optional))
+              tempVector[0] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &tempVector[1], 1,optional))
+              tempVector[1] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &tempVector[2], 1,optional))
+              tempVector[2] = 0.;
+            neues->setAtomicVelocity(repetition, tempVector);
 
             // parse forces if present
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
-              value[0] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
-              value[1] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
-              value[2] = 0.;
-            for (int d=0;d<NDIM;d++)
-              neues->Trajectory.F.at(repetition)[d] = value[d];
+            tempVector.Zero();
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &tempVector[0], 1,optional))
+              tempVector[0] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &tempVector[1], 1,optional))
+              tempVector[1] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &tempVector[2], 1,optional))
+              tempVector[2] = 0.;
+            neues->setAtomicForce(repetition, tempVector);
 
   //            Log() << Verbose(0) << "Parsed position of step " << (repetition) << ": (";
-  //            for (int d=0;d<NDIM;d++)
-  //              Log() << Verbose(0) << neues->Trajectory.R.at(repetition).x[d] << " ";          // next step
-  //            Log() << Verbose(0) << ")\t(";
-  //            for (int d=0;d<NDIM;d++)
-  //              Log() << Verbose(0) << neues->Trajectory.U.at(repetition).x[d] << " ";          // next step
-  //            Log() << Verbose(0) << ")\t(";
-  //            for (int d=0;d<NDIM;d++)
-  //              Log() << Verbose(0) << neues->Trajectory.F.at(repetition).x[d] << " ";          // next step
+  //            Log() << Verbose(0) << neues->getPosition(repetition) << " ";          // next step
+  //            Log() << Verbose(0) << neues->getAtomicVelocity(repetition) << " ";          // next step
+  //            Log() << Verbose(0) << neues->getAtomicForce(repetition) << " ";          // next step
   //            Log() << Verbose(0) << ")" << endl;
           }
@@ -603,15 +597,21 @@
             neues = AtomList[i][j];
           // then parse for each atom the coordinates as often as present
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &position[0], repetition,critical);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], repetition,critical);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], repetition,critical);
-          neues->setPosition(position);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->getAtomicVelocity()[0], repetition,optional))
-            neues->getAtomicVelocity()[0] = 0.;
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->getAtomicVelocity()[1], repetition,optional))
-            neues->getAtomicVelocity()[1] = 0.;
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->getAtomicVelocity()[2], repetition,optional))
-            neues->getAtomicVelocity()[2] = 0.;
+          tempVector.Zero();
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &tempVector[0], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &tempVector[1], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &tempVector[2], repetition,critical);
+          neues->setPosition(tempVector);
+          // fixed ion
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &_fixedion, repetition,critical);
+          neues->setFixedIon(_fixedion);
+          // and velocity
+          tempVector.Zero();
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &tempVector[0], repetition,optional))
+            tempVector[0] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &tempVector[1], repetition,optional))
+            tempVector[1] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &tempVector[2], repetition,optional))
+            tempVector[2] = 0.;
+          neues->setAtomicVelocity(tempVector);
           // here we don't care if forces are present (last in trajectories is always equal to current position)
           neues->setType(elementhash[i]); // find element type
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/molecule.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -365,5 +365,5 @@
       FirstOtherAtom->setType(1);  // element is Hydrogen
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
+      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
@@ -430,7 +430,7 @@
       SecondOtherAtom->setType(1);
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
+      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
+      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -485,9 +485,9 @@
       ThirdOtherAtom->setType(1);
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
+      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
+      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
+      ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -582,5 +582,5 @@
   istringstream *input = NULL;
   int NumberOfAtoms = 0; // atom number in xyz read
-  int i, j; // loop variables
+  int i; // loop variables
   atom *Walker = NULL;  // pointer to added atom
   char shorthand[3];  // shorthand for atom name
@@ -618,15 +618,11 @@
       Walker->setType(1);
     }
-    if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
-      Walker->Trajectory.R.resize(MDSteps+10);
-      Walker->Trajectory.U.resize(MDSteps+10);
-      Walker->Trajectory.F.resize(MDSteps+10);
+    if (Walker->getTrajectorySize() <= (unsigned int)MDSteps) {
+      Walker->ResizeTrajectory(MDSteps+10);
     }
     Walker->setPosition(Vector(x));
-    for(j=NDIM;j--;) {
-      Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
-      Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
-      Walker->Trajectory.F.at(MDSteps-1)[j] = 0;
-    }
+    Walker->setPosition(MDSteps-1, Vector(x));
+    Walker->setAtomicVelocity(MDSteps-1, zeroVec);
+    Walker->setAtomicForce(MDSteps-1, zeroVec);
     AddAtom(Walker);  // add to molecule
     delete(item);
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/molecule_dynamics.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -60,17 +60,17 @@
     // determine normalized trajectories direction vector (n1, n2)
     Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
-    trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);
+    trajectory1 = Sprinter->getPosition(Params.endstep) - Walker->getPosition(Params.startstep);
     trajectory1.Normalize();
     Norm1 = trajectory1.Norm();
     Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
-    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
+    trajectory2 = Sprinter->getPosition(Params.endstep) - (*iter)->getPosition(Params.startstep);
     trajectory2.Normalize();
     Norm2 = trajectory1.Norm();
     // check whether either is zero()
     if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
-      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
+      tmp = Walker->getPosition(Params.startstep).distance((*iter)->getPosition(Params.startstep));
     } else if (Norm1 < MYEPSILON) {
       Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
-      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
+      trajectory1 = Sprinter->getPosition(Params.endstep) - (*iter)->getPosition(Params.startstep);
       trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
       trajectory1 -= trajectory2;   // project the part in norm direction away
@@ -78,5 +78,5 @@
     } else if (Norm2 < MYEPSILON) {
       Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
-      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
+      trajectory2 = Sprinter->getPosition(Params.endstep) - Walker->getPosition(Params.startstep);  // copy second offset
       trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
       trajectory2 -= trajectory1;   // project the part in norm direction away
@@ -87,5 +87,5 @@
   //        Log() << Verbose(0) << " and ";
   //        Log() << Verbose(0) << trajectory2;
-      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
+      tmp = Walker->getPosition(Params.startstep).distance((*iter)->getPosition(Params.startstep));
   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
     } else { // determine distance by finding minimum distance
@@ -106,5 +106,5 @@
         gsl_matrix_set(A, 1, i, trajectory2[i]);
         gsl_matrix_set(A, 2, i, normal[i]);
-        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
+        gsl_vector_set(x,i, (Walker->getPosition(Params.startstep)[i] - (*iter)->getPosition(Params.startstep)[i]));
       }
       // solve the linear system by Householder transformations
@@ -117,6 +117,6 @@
       trajectory2.Scale(gsl_vector_get(x,1));
       normal.Scale(gsl_vector_get(x,2));
-      TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
-                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
+      TestVector = (*iter)->getPosition(Params.startstep) + trajectory2 + normal
+                   - (Walker->getPosition(Params.startstep) + trajectory1);
       if (TestVector.Norm() < MYEPSILON) {
   //          Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
@@ -150,5 +150,5 @@
   //    atom *Sprinter = PermutationMap[Walker->nr];
   //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
-  //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
+  //        Log() << Verbose(0) << Sprinter->getPosition(endstep);
   //        Log() << Verbose(0) << ", penalting." << endl;
       result += Params.PenaltyConstants[2];
@@ -183,5 +183,5 @@
     // first term: distance to target
     Runner = Params.PermutationMap[(*iter)->nr];   // find target point
-    tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
+    tmp = ((*iter)->getPosition(Params.startstep).distance(Runner->getPosition(Params.endstep)));
     tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     result += Params.PenaltyConstants[0] * tmp;
@@ -239,5 +239,5 @@
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
-      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).distance((*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
+      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->getPosition(Params.startstep).distance((*runner)->getPosition(Params.endstep)), (*runner)) );
     }
   }
@@ -489,5 +489,5 @@
     // set forces
     for (int i=NDIM;i++;)
-      Force->Matrix[0][_atom->nr][5+i] += 2.*constant*sqrt(_atom->Trajectory.R.at(startstep).distance(Sprinter->Trajectory.R.at(endstep)));
+      Force->Matrix[0][_atom->nr][5+i] += 2.*constant*sqrt(_atom->getPosition(startstep).distance(Sprinter->getPosition(endstep)));
   }
   DoLog(1) && (Log() << Verbose(1) << "done." << endl);
@@ -525,5 +525,6 @@
 
   // check whether we have sufficient space in Trajectories for each atom
-  for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps));
+  LOG(1, "STATUS: Extending each trajectory size to " << MaxSteps+1 << ".");
+  for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps+1));
   // push endstep to last one
   for_each(atoms.begin(),atoms.end(),boost::bind(&atom::CopyStepOnStep,_1,MaxSteps,endstep));
@@ -538,14 +539,10 @@
       // add to molecule list
       Sprinter = mol->AddCopyAtom((*iter));
-      for (int n=NDIM;n--;) {
-        Sprinter->set(n, (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps));
-        // add to Trajectories
-        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
-        if (step < MaxSteps) {
-          (*iter)->Trajectory.R.at(step)[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
-          (*iter)->Trajectory.U.at(step)[n] = 0.;
-          (*iter)->Trajectory.F.at(step)[n] = 0.;
-        }
-      }
+      // add to Trajectories
+      Vector temp = (*iter)->getPosition(startstep) + (PermutationMap[(*iter)->nr]->getPosition(endstep) - (*iter)->getPosition(startstep))*((double)step/(double)MaxSteps);
+      Sprinter->setPosition(temp);
+      (*iter)->setAtomicVelocity(step, zeroVec);
+      (*iter)->setAtomicForce(step, zeroVec);
+      //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
     }
   }
@@ -620,7 +617,7 @@
   // and perform Verlet integration for each atom with position, velocity and force vector
   // check size of vectors
-  for_each(atoms.begin(),
-           atoms.end(),
-           boost::bind(&atom::VelocityVerletUpdate,_1,MDSteps+1, &configuration, &Force, (const size_t) 0));
+  BOOST_FOREACH(atom *_atom, atoms) {
+    _atom->VelocityVerletUpdate(_atom->nr, MDSteps+1, &configuration, &Force, (const size_t) 0);
+  }
 
   // correct velocities (rather momenta) so that center of mass remains motionless
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision 54b42e9e210027bb6e47d8594317475629bdb453)
+++ src/molecule_geometry.cpp	(revision 6625c3d80c1d7b0bad97e9d543982f10829d358b)
@@ -58,6 +58,8 @@
   // go through all atoms
   BOOST_FOREACH(atom* iter, atoms){
+    std::cout << "atom before is at " << *iter << std::endl;
     *iter -= *Center;
-    *iter -= *CenterBox;
+    *iter += *CenterBox;
+    std::cout << "atom after is at " << *iter << std::endl;
   }
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
@@ -314,7 +316,9 @@
 {
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    for (int j=0;j<MDSteps;j++)
-      (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
-    (*iter)->ScaleAll(*factor);
+    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
+      Vector temp = (*iter)->getPosition(j);
+      temp.ScaleAll(*factor);
+      (*iter)->setPosition(j,temp);
+    }
   }
 };
@@ -326,7 +330,9 @@
 {
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    for (int j=0;j<MDSteps;j++)
-      (*iter)->Trajectory.R.at(j) += (*trans);
-    *(*iter) += (*trans);
+    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
+      Vector temp = (*iter)->getPosition(j);
+      temp += (*trans);
+      (*iter)->setPosition(j, temp);
+    }
   }
 };
@@ -334,5 +340,5 @@
 /** Translate the molecule periodically in the box.
  * \param trans[] translation vector.
- * TODO treatment of trajetories missing
+ * TODO treatment of trajectories missing
  */
 void molecule::TranslatePeriodically(const Vector *trans)
@@ -438,7 +444,8 @@
     (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
     for (int j=0;j<MDSteps;j++) {
-      tmp = (*iter)->Trajectory.R.at(j)[0];
-      (*iter)->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
-      (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
+      Vector temp;
+      temp[0] =  cos(alpha) * (*iter)->getPosition(j)[0] + sin(alpha) * (*iter)->getPosition(j)[2];
+      temp[2] = -sin(alpha) * (*iter)->getPosition(j)[0] + cos(alpha) * (*iter)->getPosition(j)[2];
+      (*iter)->setPosition(j,temp);
     }
   }
@@ -457,7 +464,8 @@
     (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
     for (int j=0;j<MDSteps;j++) {
-      tmp = (*iter)->Trajectory.R.at(j)[1];
-      (*iter)->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
-      (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
+      Vector temp;
+      temp[1] =  cos(alpha) * (*iter)->getPosition(j)[1] + sin(alpha) * (*iter)->getPosition(j)[2];
+      temp[2] = -sin(alpha) * (*iter)->getPosition(j)[1] + cos(alpha) * (*iter)->getPosition(j)[2];
+      (*iter)->setPosition(j,temp);
     }
   }
