Index: molecuilder/src/atom.cpp
===================================================================
--- molecuilder/src/atom.cpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/atom.cpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -64,4 +64,7 @@
   Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
   Free<char>(&Name, "atom::~atom: *Name");
+  Trajectory.R.clear();
+  Trajectory.U.clear();
+  Trajectory.F.clear();
 };
 
@@ -118,5 +121,5 @@
  * \param *comment commentary after '#' sign
  */
-bool atom::Output(int ElementNo, int AtomNo, ofstream *out, const char *comment) const
+bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const
 {
   if (out != NULL) {
@@ -134,5 +137,5 @@
     return false;
 };
-bool atom::Output(int *ElementNo, int *AtomNo, ofstream *out, const char *comment)
+bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
 {
   AtomNo[type->Z]++;  // increment number
@@ -164,4 +167,24 @@
 };
 
+/** Output of a single atom as one lin in xyz file.
+ * \param *out stream to output to
+ */
+bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const
+{
+  AtomNo[type->Z]++;
+  if (out != NULL) {
+    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
+    *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
+    *out << "\t" << FixedIon;
+    if (Trajectory.U.at(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t";
+    if (Trajectory.F.at(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[2] << "\t";
+    *out << "\t# Number in molecule " << nr << endl;
+    return true;
+  } else
+    return false;
+};
+
 ostream & operator << (ostream &ost, const atom &a)
 {
Index: molecuilder/src/atom.hpp
===================================================================
--- molecuilder/src/atom.hpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/atom.hpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -18,4 +18,5 @@
 
 #include <iostream>
+#include <vector>
 
 #include "element.hpp"
@@ -28,6 +29,14 @@
 class atom : public TesselPoint {
   public:
-    Vector x;       //!< coordinate array of atom, giving position within cell
-    Vector v;       //!< velocity array of atom
+    struct
+    {
+      vector<Vector> R;  //!< position vector
+      vector<Vector> U;  //!< velocity vector
+      vector<Vector> F;  //!< last force vector
+    } Trajectory;
+
+    Vector x;       //!< coordinate vector of atom, giving last position within cell
+    Vector v;       //!< velocity vector of atom, giving last velocity within cell
+    Vector F;       //!< Force vector of atom, giving last force within cell
     element *type;  //!< pointing to element
     atom *previous; //!< previous atom in molecule list
@@ -51,7 +60,8 @@
   virtual ~atom();
 
-  bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
-  bool Output(int *ElementNo, int *AtomNo, ofstream *out, const char *comment = NULL);
+  bool Output(ofstream *out, int ElementNo, int AtomNo, const char *comment = NULL) const;
+  bool Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment = NULL);
   bool OutputXYZLine(ofstream *out) const;
+  bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
   void EqualsFather ( atom *ptr, atom **res );
   void CorrectFather();
@@ -66,4 +76,5 @@
 };
 
+
 ostream & operator << (ostream &ost, const atom &a);
 
Index: molecuilder/src/boundary.cpp
===================================================================
--- molecuilder/src/boundary.cpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/boundary.cpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -1029,5 +1029,5 @@
 
   // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
-  mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
+  //mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
 //  mol->GoToFirst();
 //  class TesselPoint *Runner = NULL;
Index: molecuilder/src/boundary.hpp
===================================================================
--- molecuilder/src/boundary.hpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/boundary.hpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -16,5 +16,5 @@
 
 #define DEBUG 1
-#define DoSingleStepOutput 1
+#define DoSingleStepOutput 0
 #define SingleStepWidth 1
 
Index: molecuilder/src/config.cpp
===================================================================
--- molecuilder/src/config.cpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/config.cpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -964,14 +964,14 @@
 
             // check size of vectors
-            if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
+            if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
               //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
-              mol->Trajectories[neues].R.resize(repetition+10);
-              mol->Trajectories[neues].U.resize(repetition+10);
-              mol->Trajectories[neues].F.resize(repetition+10);
+              neues->Trajectory.R.resize(repetition+10);
+              neues->Trajectory.U.resize(repetition+10);
+              neues->Trajectory.F.resize(repetition+10);
             }
 
             // put into trajectories list
             for (int d=0;d<NDIM;d++)
-              mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
+              neues->Trajectory.R.at(repetition).x[d] = neues->x.x[d];
 
             // parse velocities if present
@@ -983,5 +983,5 @@
               neues->v.x[2] = 0.;
             for (int d=0;d<NDIM;d++)
-              mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
+              neues->Trajectory.U.at(repetition).x[d] = neues->v.x[d];
 
             // parse forces if present
@@ -993,15 +993,15 @@
               value[2] = 0.;
             for (int d=0;d<NDIM;d++)
-              mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
+              neues->Trajectory.F.at(repetition).x[d] = value[d];
 
   //            cout << "Parsed position of step " << (repetition) << ": (";
   //            for (int d=0;d<NDIM;d++)
-  //              cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";          // next step
+  //              cout << neues->Trajectory.R.at(repetition).x[d] << " ";          // next step
   //            cout << ")\t(";
   //            for (int d=0;d<NDIM;d++)
-  //              cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";          // next step
+  //              cout << neues->Trajectory.U.at(repetition).x[d] << " ";          // next step
   //            cout << ")\t(";
   //            for (int d=0;d<NDIM;d++)
-  //              cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";          // next step
+  //              cout << neues->Trajectory.F.at(repetition).x[d] << " ";          // next step
   //            cout << ")" << endl;
           }
Index: molecuilder/src/molecule.cpp
===================================================================
--- molecuilder/src/molecule.cpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/molecule.cpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -458,14 +458,14 @@
       Walker->type = elemente->FindElement(1);
     }
-    if (Trajectories[Walker].R.size() <= (unsigned int)MDSteps) {
-      Trajectories[Walker].R.resize(MDSteps+10);
-      Trajectories[Walker].U.resize(MDSteps+10);
-      Trajectories[Walker].F.resize(MDSteps+10);
+    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);
     }
     for(j=NDIM;j--;) {
       Walker->x.x[j] = x[j];
-      Trajectories[Walker].R.at(MDSteps-1).x[j] = x[j];
-      Trajectories[Walker].U.at(MDSteps-1).x[j] = 0;
-      Trajectories[Walker].F.at(MDSteps-1).x[j] = 0;
+      Walker->Trajectory.R.at(MDSteps-1).x[j] = x[j];
+      Walker->Trajectory.U.at(MDSteps-1).x[j] = 0;
+      Walker->Trajectory.F.at(MDSteps-1).x[j] = 0;
     }
     AddAtom(Walker);  // add to molecule
@@ -624,5 +624,4 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
-  Trajectories.erase(pointer);
   return remove(pointer, start, end);
 };
@@ -642,5 +641,4 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
-  Trajectories.erase(pointer);
   unlink(pointer);
   return true;
@@ -726,5 +724,5 @@
         ElementNo[i] = current++;
     }
-    ActOnAllAtoms( &atom::Output, ElementNo, AtomNo, out, (const char *) NULL ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
+    ActOnAllAtoms( &atom::Output, out, ElementNo, AtomNo, (const char *) NULL ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
     return true;
   }
@@ -759,17 +757,5 @@
           ElementNo[i] = current++;
       }
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        AtomNo[walker->type->Z]++;
-        *out << "Ion_Type" << ElementNo[walker->type->Z] << "_" << AtomNo[walker->type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
-        *out << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2];
-        *out << "\t" << walker->FixedIon;
-        if (Trajectories[walker].U.at(step).Norm() > MYEPSILON)
-          *out << "\t" << scientific << setprecision(6) << Trajectories[walker].U.at(step).x[0] << "\t" << Trajectories[walker].U.at(step).x[1] << "\t" << Trajectories[walker].U.at(step).x[2] << "\t";
-        if (Trajectories[walker].F.at(step).Norm() > MYEPSILON)
-          *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t";
-        *out << "\t# Number in molecule " << walker->nr << endl;
-      }
+      ActOnAllAtoms( &atom::OutputTrajectory, out, ElementNo, AtomNo, step ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
     }
     return true;
@@ -829,5 +815,5 @@
       while (walker->next != end) { // go through every atom of this element
         walker = walker->next;
-        *out << walker->type->symbol << "\t" << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2] << endl;
+        *out << walker->type->symbol << "\t" << walker->Trajectory.R.at(step).x[0] << "\t" << walker->Trajectory.R.at(step).x[1] << "\t" << walker->Trajectory.R.at(step).x[2] << endl;
       }
     }
@@ -1237,5 +1223,5 @@
       Walker = Walker->next;
       for (int i=NDIM;i--;)
-        temperature += Walker->type->mass * Trajectories[Walker].U.at(step).x[i]* Trajectories[Walker].U.at(step).x[i];
+        temperature += Walker->type->mass * Walker->Trajectory.U.at(step).x[i]* Walker->Trajectory.U.at(step).x[i];
     }
     *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
Index: molecuilder/src/molecule.hpp
===================================================================
--- molecuilder/src/molecule.hpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/molecule.hpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -65,12 +65,4 @@
 {
   bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
-};
-
-struct Trajectory
-{
-  vector<Vector> R;  //!< position vector
-  vector<Vector> U;  //!< velocity vector
-  vector<Vector> F;  //!< last force vector
-  atom *ptr;         //!< pointer to atom whose trajectory we contain
 };
 
@@ -101,5 +93,4 @@
     bond *last;         //!< end of bond list
     bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
-    map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
     int MDSteps;        //!< The number of MD steps in Trajectories
     int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
Index: molecuilder/src/molecule_dynamics.cpp
===================================================================
--- molecuilder/src/molecule_dynamics.cpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/molecule_dynamics.cpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -46,5 +46,5 @@
     // first term: distance to target
     Runner = PermutationMap[Walker->nr];   // find target point
-    tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
+    tmp = (Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(endstep)));
     tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     result += constants[0] * tmp;
@@ -59,20 +59,20 @@
       // determine normalized trajectories direction vector (n1, n2)
       Sprinter = PermutationMap[Walker->nr];   // find first target point
-      trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
-      trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
+      trajectory1.CopyVector(&Sprinter->Trajectory.R.at(endstep));
+      trajectory1.SubtractVector(&Walker->Trajectory.R.at(startstep));
       trajectory1.Normalize();
       Norm1 = trajectory1.Norm();
       Sprinter = PermutationMap[Runner->nr];   // find second target point
-      trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
-      trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
+      trajectory2.CopyVector(&Sprinter->Trajectory.R.at(endstep));
+      trajectory2.SubtractVector(&Runner->Trajectory.R.at(startstep));
       trajectory2.Normalize();
       Norm2 = trajectory1.Norm();
       // check whether either is zero()
       if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
-        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
+        tmp = Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(startstep));
       } else if (Norm1 < MYEPSILON) {
         Sprinter = PermutationMap[Walker->nr];   // find first target point
-        trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
-        trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
+        trajectory1.CopyVector(&Sprinter->Trajectory.R.at(endstep));  // copy first offset
+        trajectory1.SubtractVector(&Runner->Trajectory.R.at(startstep));  // subtract second offset
         trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
         trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
@@ -80,6 +80,6 @@
       } else if (Norm2 < MYEPSILON) {
         Sprinter = PermutationMap[Runner->nr];   // find second target point
-        trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
-        trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
+        trajectory2.CopyVector(&Sprinter->Trajectory.R.at(endstep));  // copy second offset
+        trajectory2.SubtractVector(&Walker->Trajectory.R.at(startstep));  // subtract first offset
         trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
         trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
@@ -90,5 +90,5 @@
 //        *out << " and ";
 //        *out << trajectory2;
-        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
+        tmp = Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(startstep));
 //        *out << " with distance " << tmp << "." << endl;
       } else { // determine distance by finding minimum distance
@@ -109,5 +109,5 @@
           gsl_matrix_set(A, 1, i, trajectory2.x[i]);
           gsl_matrix_set(A, 2, i, normal.x[i]);
-          gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
+          gsl_vector_set(x,i, (Walker->Trajectory.R.at(startstep).x[i] - Runner->Trajectory.R.at(startstep).x[i]));
         }
         // solve the linear system by Householder transformations
@@ -117,10 +117,10 @@
 //        *out << " with distance " << tmp << "." << endl;
         // test whether we really have the intersection (by checking on c_1 and c_2)
-        TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
+        TestVector.CopyVector(&Runner->Trajectory.R.at(startstep));
         trajectory2.Scale(gsl_vector_get(x,1));
         TestVector.AddVector(&trajectory2);
         normal.Scale(gsl_vector_get(x,2));
         TestVector.AddVector(&normal);
-        TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
+        TestVector.SubtractVector(&Walker->Trajectory.R.at(startstep));
         trajectory1.Scale(gsl_vector_get(x,0));
         TestVector.SubtractVector(&trajectory1);
@@ -148,5 +148,5 @@
         Sprinter = PermutationMap[Walker->nr];
 //        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
-//        *out << Trajectories[Sprinter].R.at(endstep);
+//        *out << Sprinter->Trajectory.R.at(endstep);
 //        *out << ", penalting." << endl;
         result += constants[2];
@@ -240,5 +240,5 @@
     while(Runner->next != end) {
       Runner = Runner->next;
-      DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
+      DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(endstep)), Runner) );
     }
   }
@@ -398,5 +398,5 @@
     // set forces
     for (int i=NDIM;i++;)
-      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
+      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Walker->Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
   }
   *out << Verbose(1) << "done." << endl;
@@ -436,9 +436,9 @@
   while (Walker->next != end) {
     Walker = Walker->next;
-    if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
+    if (Walker->Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
       //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
-      Trajectories[Walker].R.resize(MaxSteps+1);
-      Trajectories[Walker].U.resize(MaxSteps+1);
-      Trajectories[Walker].F.resize(MaxSteps+1);
+      Walker->Trajectory.R.resize(MaxSteps+1);
+      Walker->Trajectory.U.resize(MaxSteps+1);
+      Walker->Trajectory.F.resize(MaxSteps+1);
     }
   }
@@ -448,7 +448,7 @@
     Walker = Walker->next;
     for (int n=NDIM;n--;) {
-      Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
-      Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
-      Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
+      Walker->Trajectory.R.at(MaxSteps).x[n] = Walker->Trajectory.R.at(endstep).x[n];
+      Walker->Trajectory.U.at(MaxSteps).x[n] = Walker->Trajectory.U.at(endstep).x[n];
+      Walker->Trajectory.F.at(MaxSteps).x[n] = Walker->Trajectory.F.at(endstep).x[n];
     }
   }
@@ -466,11 +466,11 @@
       Sprinter = mol->AddCopyAtom(Walker);
       for (int n=NDIM;n--;) {
-        Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+        Sprinter->x.x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
         // add to Trajectories
         //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
         if (step < MaxSteps) {
-          Trajectories[Walker].R.at(step).x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
-          Trajectories[Walker].U.at(step).x[n] = 0.;
-          Trajectories[Walker].F.at(step).x[n] = 0.;
+          Walker->Trajectory.R.at(step).x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+          Walker->Trajectory.U.at(step).x[n] = 0.;
+          Walker->Trajectory.F.at(step).x[n] = 0.;
         }
       }
@@ -554,29 +554,29 @@
       //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
       // check size of vectors
-      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
+      if (walker->Trajectory.R.size() <= (unsigned int)(MDSteps)) {
         //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-        Trajectories[walker].R.resize(MDSteps+10);
-        Trajectories[walker].U.resize(MDSteps+10);
-        Trajectories[walker].F.resize(MDSteps+10);
+        walker->Trajectory.R.resize(MDSteps+10);
+        walker->Trajectory.U.resize(MDSteps+10);
+        walker->Trajectory.F.resize(MDSteps+10);
       }
 
       // Update R (and F)
       for (int d=0; d<NDIM; d++) {
-        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
-        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-        Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+        walker->Trajectory.F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+        walker->Trajectory.R.at(MDSteps).x[d] = walker->Trajectory.R.at(MDSteps-1).x[d];
+        walker->Trajectory.R.at(MDSteps).x[d] += configuration.Deltat*(walker->Trajectory.U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+        walker->Trajectory.R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(walker->Trajectory.F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
       }
       // Update U
       for (int d=0; d<NDIM; d++) {
-        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
-        Trajectories[walker].U.at(MDSteps).x[d] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
+        walker->Trajectory.U.at(MDSteps).x[d] = walker->Trajectory.U.at(MDSteps-1).x[d];
+        walker->Trajectory.U.at(MDSteps).x[d] += configuration.Deltat * (walker->Trajectory.F.at(MDSteps).x[d]+walker->Trajectory.F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
       }
 //      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
 //      for (int d=0;d<NDIM;d++)
-//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
+//        out << walker->Trajectory.R.at(MDSteps).x[d] << " ";          // next step
 //      out << ")\t(";
 //      for (int d=0;d<NDIM;d++)
-//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
+//        cout << walker->Trajectory.U.at(MDSteps).x[d] << " ";          // next step
 //      out << ")" << endl;
             // next atom
@@ -592,5 +592,5 @@
     IonMass += walker->type->mass;  // sum up total mass
     for(int d=0;d<NDIM;d++) {
-      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
+      Vector[d] += walker->Trajectory.U.at(MDSteps).x[d]*walker->type->mass;
     }
   }
@@ -603,6 +603,6 @@
     walker = walker->next;
     for(int d=0;d<NDIM;d++) {
-      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
-      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
+      walker->Trajectory.U.at(MDSteps).x[d] -= Vector[d];
+      ActualTemp += 0.5 * walker->type->mass * walker->Trajectory.U.at(MDSteps).x[d] * walker->Trajectory.U.at(MDSteps).x[d];
     }
   }
@@ -663,5 +663,5 @@
           walker = walker->next;
           IonMass = walker->type->mass;
-          U = Trajectories[walker].U.at(MDSteps).x;
+          U = walker->Trajectory.U.at(MDSteps).x;
           if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
             for (d=0; d<NDIM; d++) {
@@ -678,6 +678,6 @@
         walker = walker->next;
         IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
+        U = walker->Trajectory.U.at(MDSteps).x;
+        F = walker->Trajectory.F.at(MDSteps).x;
         if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
           for (d=0; d<NDIM; d++) {
@@ -691,6 +691,6 @@
         walker = walker->next;
         IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
+        U = walker->Trajectory.U.at(MDSteps).x;
+        F = walker->Trajectory.F.at(MDSteps).x;
         if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
           for (d=0; d<NDIM; d++) {
@@ -713,6 +713,6 @@
         IonMass = walker->type->mass;
         sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
+        U = walker->Trajectory.U.at(MDSteps).x;
+        F = walker->Trajectory.F.at(MDSteps).x;
         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
           // throw a dice to determine whether it gets hit by a heat bath particle
@@ -736,6 +736,6 @@
         walker = walker->next;
         IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
+        U = walker->Trajectory.U.at(MDSteps).x;
+        F = walker->Trajectory.F.at(MDSteps).x;
         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
           for (d=0; d<NDIM; d++) {
@@ -754,5 +754,5 @@
         walker = walker->next;
         IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
+        U = walker->Trajectory.U.at(MDSteps).x;
         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
           for (d=0; d<NDIM; d++) {
@@ -769,5 +769,5 @@
         walker = walker->next;
         IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
+        U = walker->Trajectory.U.at(MDSteps).x;
         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
           for (d=0; d<NDIM; d++) {
Index: molecuilder/src/molecule_geometry.cpp
===================================================================
--- molecuilder/src/molecule_geometry.cpp	(revision 8ffe326c908c4250db92da4fac0102a8a5948f68)
+++ molecuilder/src/molecule_geometry.cpp	(revision 567b7f4e79b7999319a380200ae0622682af97ec)
@@ -200,5 +200,5 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      Trajectories[ptr].R.at(j).Scale(factor);
+      ptr->Trajectory.R.at(j).Scale(factor);
     ptr->x.Scale(factor);
   }
@@ -215,5 +215,5 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      Trajectories[ptr].R.at(j).Translate(trans);
+      ptr->Trajectory.R.at(j).Translate(trans);
     ptr->x.Translate(trans);
   }
@@ -438,7 +438,7 @@
     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     for (int j=0;j<MDSteps;j++) {
-      tmp = Trajectories[ptr].R.at(j).x[0];
-      Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
-      Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
+      tmp = ptr->Trajectory.R.at(j).x[0];
+      ptr->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
+      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
     }
   }
@@ -461,7 +461,7 @@
     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     for (int j=0;j<MDSteps;j++) {
-      tmp = Trajectories[ptr].R.at(j).x[1];
-      Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
-      Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
+      tmp = ptr->Trajectory.R.at(j).x[1];
+      ptr->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
+      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
     }
   }
