Index: src/Atom/atom_atominfo.cpp
===================================================================
--- src/Atom/atom_atominfo.cpp	(revision 7e51e12e82c5637d2660d794f1e024fda52a5aa6)
+++ src/Atom/atom_atominfo.cpp	(revision b119a092aa1848936f8c80974f857ded3c098885)
@@ -3,4 +3,5 @@
  * Description: creates and alters molecular systems
  * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
+ * Copyright (C)  2014 Frederik Heber. All rights reserved.
  * 
  *
@@ -55,12 +56,8 @@
   charge(0.)
 {
-  AtomicPosition.reserve(1);
-  AtomicPosition.push_back(zeroVec);
-  AtomicVelocity.reserve(1);
-  AtomicVelocity.push_back(zeroVec);
-  AtomicForce.reserve(1);
-  AtomicForce.push_back(zeroVec);
-
-};
+  AtomicPosition.insert( std::make_pair(0, zeroVec) );
+  AtomicVelocity.insert( std::make_pair(0, zeroVec) );
+  AtomicForce.insert( std::make_pair(0, zeroVec) );
+}
 
 /** Copy constructor of class AtomInfo.
@@ -68,13 +65,11 @@
 AtomInfo::AtomInfo(const AtomInfo &_atom) :
     AtomicPosition(_atom.AtomicPosition),
+    AtomicVelocity(_atom.AtomicVelocity),
+    AtomicForce(_atom.AtomicForce),
     AtomicElement(_atom.AtomicElement),
     FixedIon(_atom.FixedIon),
     charge(_atom.charge)
 {
-  AtomicVelocity.reserve(1);
-  AtomicVelocity.push_back(zeroVec);
-  AtomicForce.reserve(1);
-  AtomicForce.push_back(zeroVec);
-};
+}
 
 AtomInfo::AtomInfo(const VectorInterface &_v) :
@@ -83,11 +78,7 @@
     charge(0.)
 {
-  if (AtomicPosition.size() < 1)
-    AtomicPosition.resize(1);
-  AtomicPosition[0] = _v.getPosition();
-  AtomicVelocity.reserve(1);
-  AtomicVelocity.push_back(zeroVec);
-  AtomicForce.reserve(1);
-  AtomicForce.push_back(zeroVec);
+  AtomicPosition.insert( std::make_pair(0, _v.getPosition()) );
+  AtomicVelocity.insert( std::make_pair(0, zeroVec) );
+  AtomicForce.insert( std::make_pair(0, zeroVec) );
 };
 
@@ -98,10 +89,11 @@
 };
 
-void AtomInfo::AppendTrajectoryStep()
-{
-  NOTIFY(TrajectoryChanged);
-  AtomicPosition.push_back(zeroVec);
-  AtomicVelocity.push_back(zeroVec);
-  AtomicForce.push_back(zeroVec);
+void AtomInfo::AppendTrajectoryStep(const unsigned int _step)
+{
+  if (WorldTime::getTime() == _step)
+    NOTIFY(TrajectoryChanged);
+  AtomicPosition.insert( std::make_pair(_step, zeroVec) );
+  AtomicVelocity.insert( std::make_pair(_step, zeroVec) );
+  AtomicForce.insert( std::make_pair(_step, zeroVec) );
   LOG(5,"AtomInfo::AppendTrajectoryStep() called, size is ("
       << AtomicPosition.size() << ","
@@ -110,10 +102,11 @@
 }
 
-void AtomInfo::removeTrajectoryStep()
-{
-  NOTIFY(TrajectoryChanged);
-  AtomicPosition.pop_back();
-  AtomicVelocity.pop_back();
-  AtomicForce.pop_back();
+void AtomInfo::removeTrajectoryStep(const unsigned int _step)
+{
+  if (WorldTime::getTime() == _step)
+    NOTIFY(TrajectoryChanged);
+  AtomicPosition.erase(_step);
+  AtomicVelocity.erase(_step);
+  AtomicForce.erase(_step);
   LOG(5,"AtomInfo::removeTrajectoryStep() called, size is ("
       << AtomicPosition.size() << ","
@@ -141,27 +134,19 @@
 const double& AtomInfo::operator[](size_t i) const
 {
-  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
-      "AtomInfo::operator[]() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[WorldTime::getTime()][i];
+  return atStep(i, WorldTime::getTime());
 }
 
 const double& AtomInfo::at(size_t i) const
 {
-  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
-      "AtomInfo::at() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[WorldTime::getTime()].at(i);
+  return atStep(i, WorldTime::getTime());
 }
 
 const double& AtomInfo::atStep(size_t i, unsigned int _step) const
 {
-  ASSERT(AtomicPosition.size() > _step,
-      "AtomInfo::atStep() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[_step].at(i);
+  ASSERT(!AtomicPosition.empty(),
+      "AtomInfo::operator[]() - AtomicPosition is empty.");
+  VectorTrajectory_t::const_iterator iter =
+      AtomicPosition.lower_bound(_step);
+  return iter->second[i];
 }
 
@@ -170,27 +155,45 @@
   OBSERVE;
   NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
-      "AtomInfo::set() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].at(i) = value;
+  VectorTrajectory_t::iterator iter = AtomicPosition.find(WorldTime::getTime());
+  if (iter !=  AtomicPosition.end()) {
+    iter->second[i] = value;
+  } else {
+    Vector newPos;
+    newPos[i] = value;
+#ifndef NDEBUG
+    std::pair<VectorTrajectory_t::iterator, bool> inserter =
+#endif
+        AtomicPosition.insert( std::make_pair(WorldTime::getTime(), newPos) );
+    ASSERT( inserter.second,
+        "AtomInfo::set() - time step "+toString(WorldTime::getTime())
+        +" present after all?");
+  }
+}
+
+/** Helps to determine whether the current step really exists or getPosition() has just
+ * delivered the one closest to it in the past.
+ *
+ * \param _step step to check
+ * \param true - step exists, false - step does not exist, getPosition() delivers closest
+ */
+bool AtomInfo::isStepPresent(const unsigned int _step) const
+{
+  VectorTrajectory_t::const_iterator iter =
+      AtomicPosition.find(_step);
+  return iter != AtomicPosition.end();
 }
 
 const Vector& AtomInfo::getPosition() const
 {
-  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
-      "AtomInfo::getPosition() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[WorldTime::getTime()];
+  return getPositionAtStep(WorldTime::getTime());
 }
 
 const Vector& AtomInfo::getPositionAtStep(const unsigned int _step) const
 {
-  ASSERT(_step < AtomicPosition.size(),
-      "AtomInfo::getPositionAtStep() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[_step];
+  ASSERT(!AtomicPosition.empty(),
+      "AtomInfo::operator[]() - AtomicPosition is empty.");
+  VectorTrajectory_t::const_iterator iter =
+      AtomicPosition.lower_bound(_step);
+  return iter->second;
 }
 
@@ -207,115 +210,89 @@
 {
   const element *elem = World::getInstance().getPeriode()->FindElement(Z);
-  if (elem != NULL) {
-    OBSERVE;
-    NOTIFY(AtomObservable::ElementChanged);
-    AtomicElement = Z;
-  }
-}
-
-//Vector& AtomInfo::getAtomicVelocity()
-//{
-//  return AtomicVelocity[0];
-//}
-
-//Vector& AtomInfo::getAtomicVelocity(const int _step)
-//{
-//  ASSERT(_step < AtomicVelocity.size(),
-//      "AtomInfo::getAtomicVelocity() - Access out of range.");
-//  return AtomicVelocity[_step];
-//}
+  setType(elem);
+}
 
 const Vector& AtomInfo::getAtomicVelocity() const
 {
-  ASSERT(AtomicVelocity.size() > 0,
-      "AtomInfo::getAtomicVelocity() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicVelocity[WorldTime::getTime()];
+  return getAtomicVelocityAtStep(WorldTime::getTime());
 }
 
 const Vector& AtomInfo::getAtomicVelocityAtStep(const unsigned int _step) const
 {
-  ASSERT(_step < AtomicVelocity.size(),
-      "AtomInfo::getAtomicVelocity() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicVelocity[_step];
+  ASSERT(!AtomicVelocity.empty(),
+      "AtomInfo::operator[]() - AtomicVelocity is empty.");
+  VectorTrajectory_t::const_iterator iter =
+      AtomicVelocity.lower_bound(_step);
+  // special, we only interpolate between present time steps not into the future
+  if (_step > AtomicVelocity.begin()->first)
+    return zeroVec;
+  else
+    return iter->second;
 }
 
 void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
 {
+  setAtomicVelocityAtStep(WorldTime::getTime(), _newvelocity);
+}
+
+void AtomInfo::setAtomicVelocityAtStep(const unsigned int _step, const Vector &_newvelocity)
+{
   OBSERVE;
-  NOTIFY(AtomObservable::VelocityChanged);
-  ASSERT(WorldTime::getTime() < AtomicVelocity.size(),
-      "AtomInfo::setAtomicVelocity() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicVelocity[WorldTime::getTime()] = _newvelocity;
-}
-
-void AtomInfo::setAtomicVelocityAtStep(const unsigned int _step, const Vector &_newvelocity)
-{
-  OBSERVE;
+  VectorTrajectory_t::iterator iter = AtomicVelocity.find(_step);
+  if (iter !=  AtomicVelocity.end()) {
+    iter->second = _newvelocity;
+  } else {
+#ifndef NDEBUG
+    std::pair<VectorTrajectory_t::iterator, bool> inserter =
+#endif
+        AtomicVelocity.insert( std::make_pair(_step, _newvelocity) );
+    ASSERT( inserter.second,
+        "AtomInfo::set() - time step "+toString(_step)
+        +" present after all?");
+  }
   if (WorldTime::getTime() == _step)
     NOTIFY(AtomObservable::VelocityChanged);
-  const unsigned int size = AtomicVelocity.size();
-  ASSERT(_step <= size,
-      "AtomInfo::setAtomicVelocityAtStep() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(size)+"].");
-  if(_step < size) {
-    AtomicVelocity[_step] = _newvelocity;
-  } else if (_step == size) {
-    UpdateSteps();
-    AtomicVelocity[_step] = _newvelocity;
-  }
 }
 
 const Vector& AtomInfo::getAtomicForce() const
 {
-  ASSERT(WorldTime::getTime() < AtomicForce.size(),
-      "AtomInfo::getAtomicForce() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicForce[WorldTime::getTime()];
+  return getAtomicForceAtStep(WorldTime::getTime());
 }
 
 const Vector& AtomInfo::getAtomicForceAtStep(const unsigned int _step) const
 {
-  ASSERT(_step < AtomicForce.size(),
-      "AtomInfo::getAtomicForce() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicForce[_step];
+  ASSERT(!AtomicForce.empty(),
+      "AtomInfo::operator[]() - AtomicForce is empty.");
+  VectorTrajectory_t::const_iterator iter =
+      AtomicForce.lower_bound(_step);
+  // special, we only interpolate between present time steps not into the future
+  if (_step > AtomicForce.begin()->first)
+    return zeroVec;
+  else
+    return iter->second;
 }
 
 void AtomInfo::setAtomicForce(const Vector &_newforce)
 {
+  setAtomicForceAtStep(WorldTime::getTime(), _newforce);
+}
+
+void AtomInfo::setAtomicForceAtStep(const unsigned int _step, const Vector &_newforce)
+{
   OBSERVE;
-  NOTIFY(AtomObservable::ForceChanged);
-  ASSERT(WorldTime::getTime() < AtomicForce.size(),
-      "AtomInfo::setAtomicForce() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicForce[WorldTime::getTime()] = _newforce;
-}
-
-void AtomInfo::setAtomicForceAtStep(const unsigned int _step, const Vector &_newforce)
-{
-  OBSERVE;
+  VectorTrajectory_t::iterator iter = AtomicForce.find(_step);
+  if (iter !=  AtomicForce.end()) {
+    iter->second = _newforce;
+  } else {
+#ifndef NDEBUG
+    std::pair<VectorTrajectory_t::iterator, bool> inserter =
+#endif
+        AtomicForce.insert( std::make_pair(_step, _newforce) );
+    ASSERT( inserter.second,
+        "AtomInfo::set() - time step "+toString(_step)
+        +" present after all?");
+  }
   if (WorldTime::getTime() == _step)
     NOTIFY(AtomObservable::ForceChanged);
-  const unsigned int size = AtomicForce.size();
-  ASSERT(_step <= size,
-      "AtomInfo::setAtomicForce() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(AtomicPosition.size())+"].");
-  if(_step < size) {
-    AtomicForce[_step] = _newforce;
-  } else if (_step == size) {
-    UpdateSteps();
-    AtomicForce[_step] = _newforce;
-  }
 }
 
@@ -334,42 +311,29 @@
 void AtomInfo::setPosition(const Vector& _vector)
 {
+  setPositionAtStep(WorldTime::getTime(), _vector);
+}
+
+void AtomInfo::setPositionAtStep(unsigned int _step, const Vector& _vector)
+{
   OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::setPosition() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()] = _vector;
-  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
-}
-
-void AtomInfo::setPositionAtStep(unsigned int _step, const Vector& _vector)
-{
-  OBSERVE;
+  VectorTrajectory_t::iterator iter = AtomicPosition.find(_step);
+  if (iter !=  AtomicPosition.end()) {
+    iter->second = _vector;
+  } else {
+#ifndef NDEBUG
+    std::pair<VectorTrajectory_t::iterator, bool> inserter =
+#endif
+        AtomicPosition.insert( std::make_pair(_step, _vector) );
+    ASSERT( inserter.second,
+        "AtomInfo::set() - time step "+toString(_step)
+        +" present after all?");
+  }
   if (WorldTime::getTime() == _step)
     NOTIFY(AtomObservable::PositionChanged);
-  const unsigned int size = AtomicPosition.size();
-  ASSERT(_step <= size,
-      "AtomInfo::setPosition() - Access out of range: "
-      +toString(_step)
-      +" not in [0,"+toString(size)+"].");
-  if(_step < size) {
-    AtomicPosition[_step] = _vector;
-  } else if (_step == size) {
-    UpdateSteps();
-    AtomicPosition[_step] = _vector;
-  }
-  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
 }
 
 const VectorInterface& AtomInfo::operator+=(const Vector& b)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::operator+=() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()] += b;
+  setPosition(getPosition()+b);
   return *this;
 }
@@ -377,11 +341,5 @@
 const VectorInterface& AtomInfo::operator-=(const Vector& b)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::operator-=() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()] -= b;
+  setPosition(getPosition()-b);
   return *this;
 }
@@ -389,9 +347,5 @@
 Vector const AtomInfo::operator+(const Vector& b) const
 {
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::operator+() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  Vector a(AtomicPosition[WorldTime::getTime()]);
+  Vector a(getPosition());
   a += b;
   return a;
@@ -400,9 +354,5 @@
 Vector const AtomInfo::operator-(const Vector& b) const
 {
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::operator-() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  Vector a(AtomicPosition[WorldTime::getTime()]);
+  Vector a(getPosition());
   a -= b;
   return a;
@@ -411,47 +361,25 @@
 double AtomInfo::distance(const Vector &point) const
 {
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::distance() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[WorldTime::getTime()].distance(point);
+  return getPosition().distance(point);
 }
 
 double AtomInfo::DistanceSquared(const Vector &y) const
 {
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::DistanceSquared() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return AtomicPosition[WorldTime::getTime()].DistanceSquared(y);
+  return getPosition().DistanceSquared(y);
 }
 
 double AtomInfo::distance(const VectorInterface &_atom) const
 {
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::distance() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return _atom.distance(AtomicPosition[WorldTime::getTime()]);
+  return _atom.distance(getPosition());
 }
 
 double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
 {
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::DistanceSquared() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return _atom.DistanceSquared(AtomicPosition[WorldTime::getTime()]);
+  return _atom.DistanceSquared(getPosition());
 }
 
 VectorInterface &AtomInfo::operator=(const Vector& _vector)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::operator=() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()] = _vector;
+  setPosition(_vector);
   return *this;
 }
@@ -459,66 +387,38 @@
 void AtomInfo::ScaleAll(const double *factor)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::ScaleAll() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].ScaleAll(factor);
+  Vector temp(getPosition());
+  temp.ScaleAll(factor);
+  setPosition(temp);
 }
 
 void AtomInfo::ScaleAll(const Vector &factor)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::ScaleAll() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].ScaleAll(factor);
+  Vector temp(getPosition());
+  temp.ScaleAll(factor);
+  setPosition(temp);
 }
 
 void AtomInfo::Scale(const double factor)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::Scale() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].Scale(factor);
+  Vector temp(getPosition());
+  temp.Scale(factor);
+  setPosition(temp);
 }
 
 void AtomInfo::Zero()
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::Zero() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].Zero();
+  setPosition(zeroVec);
 }
 
 void AtomInfo::One(const double one)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::One() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].One(one);
+  setPosition(Vector(one,one,one));
 }
 
 void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
 {
-  OBSERVE;
-  NOTIFY(AtomObservable::PositionChanged);
-  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
-      "AtomInfo::LinearCombinationOfVectors() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  AtomicPosition[WorldTime::getTime()].LinearCombinationOfVectors(x1,x2,x3,factors);
+  Vector newPos;
+  newPos.LinearCombinationOfVectors(x1,x2,x3,factors);
+  setPosition(newPos);
 }
 
@@ -528,33 +428,44 @@
 double AtomInfo::getKineticEnergy(const unsigned int _step) const
 {
-  ASSERT(_step < AtomicPosition.size(),
-      "AtomInfo::getKineticEnergy() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return getMass() * AtomicVelocity[_step].NormSquared();
+  return getMass() * getAtomicVelocityAtStep(_step).NormSquared();
 }
 
 Vector AtomInfo::getMomentum(const unsigned int _step) const
 {
-  ASSERT(_step < AtomicPosition.size(),
-      "AtomInfo::getMomentum() - Access out of range: "
-      +toString(WorldTime::getTime())
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  return getMass()*AtomicVelocity[_step];
-}
-
-/** Extends the trajectory STL vector to the new size.
- * Does nothing if \a MaxSteps is smaller than current size.
+  return getMass() * getAtomicVelocityAtStep(_step);
+}
+
+/** Decrease the trajectory if given \a MaxSteps is smaller.
+ * Does nothing if \a MaxSteps is larger than current size.
+ *
  * \param MaxSteps
  */
 void AtomInfo::ResizeTrajectory(size_t MaxSteps)
 {
-  for (;AtomicPosition.size() <= (unsigned int)(MaxSteps);)
-    UpdateSteps();
+  // mind the reverse ordering due to std::greater, latest time steps are at beginning
+  VectorTrajectory_t::iterator positer = AtomicPosition.lower_bound(MaxSteps);
+  if (positer != AtomicPosition.begin()) {
+    if (positer->first == MaxSteps)
+      --positer;
+    AtomicPosition.erase(AtomicPosition.begin(), positer);
+  }
+  VectorTrajectory_t::iterator veliter = AtomicVelocity.lower_bound(MaxSteps);
+  if (veliter != AtomicVelocity.begin()) {
+    if (veliter->first == MaxSteps)
+      --veliter;
+    AtomicVelocity.erase(AtomicVelocity.begin(), veliter);
+  }
+  VectorTrajectory_t::iterator forceiter = AtomicForce.lower_bound(MaxSteps);
+  if (forceiter != AtomicForce.begin()) {
+    if (forceiter->first == MaxSteps)
+      --forceiter;
+    AtomicForce.erase(AtomicForce.begin(), forceiter);
+  }
 }
 
 size_t AtomInfo::getTrajectorySize() const
 {
-  return AtomicPosition.size();
+  // mind greater comp for map here: first element is latest in time steps!
+  return AtomicPosition.begin()->first+1;
 }
 
@@ -562,4 +473,22 @@
 {
   return getType()->getMass();
+}
+
+/** Helper function to either insert or assign, depending on the element being
+ * present already.
+ *
+ * \param _trajectory vector of Vectors to assign
+ * \param dest step to insert/assign to
+ * \param _newvalue new Vector value
+ */
+void assignTrajectoryElement(
+    std::map<unsigned int, Vector, std::greater<unsigned int> > &_trajectory,
+    const unsigned int dest,
+    const Vector &_newvalue)
+{
+  std::pair<std::map<unsigned int, Vector, std::greater<unsigned int> >::iterator, bool> inserter =
+      _trajectory.insert( std::make_pair(dest, _newvalue) );
+  if (!inserter.second)
+    inserter.first->second = _newvalue;
 }
 
@@ -579,17 +508,19 @@
   }
 
-  ASSERT(dest < AtomicPosition.size(),
-      "AtomInfo::CopyStepOnStep() - destination outside of current trajectory array: "
-      +toString(dest)
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  ASSERT(src < AtomicPosition.size(),
-      "AtomInfo::CopyStepOnStep() - source outside of current trajectory array: "
-      +toString(src)
-      +" not in [0,"+toString(AtomicPosition.size())+").");
-  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];
-  }
+  VectorTrajectory_t::iterator positer = AtomicPosition.find(src);
+  ASSERT( positer != AtomicPosition.end(),
+      "AtomInfo::CopyStepOnStep() - step "
+      +toString(src)+" to copy from not present in AtomicPosition.");
+  VectorTrajectory_t::iterator veliter = AtomicVelocity.find(src);
+  ASSERT( veliter != AtomicVelocity.end(),
+      "AtomInfo::CopyStepOnStep() - step "
+      +toString(src)+" to copy from not present in AtomicVelocity.");
+  VectorTrajectory_t::iterator forceiter = AtomicForce.find(src);
+  ASSERT( forceiter != AtomicForce.end(),
+      "AtomInfo::CopyStepOnStep() - step "
+      +toString(src)+" to copy from not present in AtomicForce.");
+  assignTrajectoryElement(AtomicPosition, dest, positer->second);
+  assignTrajectoryElement(AtomicVelocity, dest, veliter->second);
+  assignTrajectoryElement(AtomicForce, dest, forceiter->second);
 };
 
@@ -656,24 +587,4 @@
 };
 
-//const AtomInfo& operator*=(AtomInfo& a, const double m)
-//{
-//  a.Scale(m);
-//  return a;
-//}
-//
-//AtomInfo const operator*(const AtomInfo& a, const double m)
-//{
-//  AtomInfo copy(a);
-//  copy *= m;
-//  return copy;
-//}
-//
-//AtomInfo const operator*(const double m, const AtomInfo& a)
-//{
-//  AtomInfo copy(a);
-//  copy *= m;
-//  return copy;
-//}
-
 std::ostream & AtomInfo::operator << (std::ostream &ost) const
 {
