Index: molecuilder/src/atom.cpp
===================================================================
--- molecuilder/src/atom.cpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/atom.cpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -16,4 +16,87 @@
 /************************************* Functions for class atom *************************************/
 
+/** Constructor of class AtomInfo.
+ */
+AtomInfo::AtomInfo() : type(NULL) {};
+
+/** Destructor of class AtomInfo.
+ */
+AtomInfo::~AtomInfo()
+{
+};
+
+/** Constructor of class TrajectoryParticleInfo.
+ */
+TrajectoryParticleInfo::TrajectoryParticleInfo() : FixedIon(0) {};
+
+/** Destructor of class TrajectoryParticleInfo.
+ */
+TrajectoryParticleInfo::~TrajectoryParticleInfo()
+{
+  Trajectory.R.clear();
+  Trajectory.U.clear();
+  Trajectory.F.clear();
+};
+
+/** Constructor of class TrajectoryParticle.
+ */
+TrajectoryParticle::TrajectoryParticle()
+{
+};
+
+/** Destructor of class TrajectoryParticle.
+ */
+TrajectoryParticle::~TrajectoryParticle()
+{
+};
+
+/** Constructor of class GraphNodeInfo.
+ */
+GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
+
+/** Destructor of class GraphNodeInfo.
+ */
+GraphNodeInfo::~GraphNodeInfo()
+{
+  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
+};
+
+/** Constructor of class GraphNode.
+ */
+GraphNode::GraphNode()
+{
+};
+
+/** Destructor of class GraphNode.
+ */
+GraphNode::~GraphNode()
+{
+};
+
+/** Constructor of class BondedParticleInfo.
+ */
+BondedParticleInfo::BondedParticleInfo() : AdaptiveOrder(0), MaxOrder(false) {};
+
+/** Destructor of class BondedParticleInfo.
+ */
+BondedParticleInfo::~BondedParticleInfo()
+{
+};
+
+/** Constructor of class BondedParticle.
+ */
+BondedParticle::BondedParticle(){};
+
+/** Destructor of class BondedParticle.
+ */
+BondedParticle::~BondedParticle()
+{
+  BondList::const_iterator Runner;
+  while (!ListOfBonds.empty()) {
+    Runner = ListOfBonds.begin();
+    removewithoutcheck(*Runner);
+  }
+};
+
 /** Constructor of class atom.
  */
@@ -23,15 +106,6 @@
   previous = NULL;
   next = NULL;
-  Ancestor = NULL;
-  type = NULL;
-  sort = NULL;
-  FixedIon = 0;
-  GraphNr = -1;
-  ComponentNr = NULL;
-  IsCyclic = false;
-  SeparationVertex = false;
-  LowpointNr = -1;
-  AdaptiveOrder = 0;
-  MaxOrder = false;
+  sort = &nr;
+
   // set LCNode::Vector to our Vector
   node = &x;
@@ -42,21 +116,12 @@
 atom::atom(atom *pointer)
 {
-  Name = NULL;
   previous = NULL;
   next = NULL;
   father = pointer;  // generally, father is itself
-  Ancestor = NULL;
-  GraphNr = -1;
-  ComponentNr = NULL;
-  IsCyclic = false;
-  SeparationVertex = false;
-  LowpointNr = -1;
-  AdaptiveOrder = 0;
-  MaxOrder = false;
+
   type = pointer->type;  // copy element of atom
   x.CopyVector(&pointer->x); // copy coordination
   v.CopyVector(&pointer->v); // copy velocity
   FixedIon = pointer->FixedIon;
-  nr = -1;
   sort = &nr;
   node = &x;
@@ -68,15 +133,5 @@
 atom::~atom()
 {
-  BondList::const_iterator Runner;
-  while (!ListOfBonds.empty()) {
-    Runner = ListOfBonds.begin();
-    removewithoutcheck(*Runner);
-  }
   unlink(this);
-  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
-  Free<char>(&Name, "atom::~atom: *Name");
-  Trajectory.R.clear();
-  Trajectory.U.clear();
-  Trajectory.F.clear();
 };
 
@@ -130,5 +185,5 @@
  * \param bonds times bond::BondDegree
  */
-int atom::CountBonds() const
+int BondedParticle::CountBonds() const
 {
   int NoBonds = 0;
@@ -233,70 +288,56 @@
 };
 
-/** Output graph info of this atom.
+/** Outputs the MPQC configuration line for this atom.
  * \param *out output stream
- */
-void atom::OutputGraphInfo(ofstream *out) const
-{
-  *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
-  OutputComponentNumber(out);
-  *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
-};
-
-/** Output a list of flags, stating whether the bond was visited or not.
- * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
- * \param *out output stream for debugging
- */
-void atom::OutputComponentNumber(ofstream *out) const
-{
-  if (ComponentNr != NULL) {
-    for (int i=0; ComponentNr[i] != -1; i++)
-      *out << ComponentNr[i] << " ";
-  }
-};
-
-
-/** Prints all bonds of this atom with total degree.
- * \param *out stream to output to
- * \return true - \a *out present, false - \a *out is NULL
- */
-bool atom::OutputBondOfAtom(ofstream *out) const
-{
-  if (out != NULL) {
-#ifdef ADDHYDROGEN
-    if (type->Z != 1) {   // regard only non-hydrogen
-#endif
-      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
-      int TotalDegree = 0;
-      for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
-        *out << **Runner << "\t";
-        TotalDegree += (*Runner)->BondDegree;
-      }
-      *out << " -- TotalDegree: " << TotalDegree << endl;
-#ifdef ADDHYDROGEN
-    }
-#endif
+ * \param *center center of molecule subtracted from position
+ * \param *AtomNo pointer to atom counter that is increased by one
+ */
+void atom::OutputMPQCLine(ofstream *out, Vector *center, int *AtomNo = NULL) const
+{
+  *out << "\t\t" << type->symbol << " [ " << x.x[0]-center->x[0] << "\t" << x.x[1]-center->x[1] << "\t" << x.x[2]-center->x[2] << " ]" << endl;
+  if (AtomNo != NULL)
+    *AtomNo++;
+};
+
+ostream & operator << (ostream &ost, const ParticleInfo &a)
+{
+  ost << "[" << a.Name << "|" << &a << "]";
+  return ost;
+};
+
+ostream & ParticleInfo::operator << (ostream &ost)
+{
+  ost << "[" << Name << "|" << this << "]";
+  return ost;
+};
+
+/** Compares the indices of \a this atom with a given \a ptr.
+ * \param ptr atom to compare index against
+ * \return true - this one's is smaller, false - not
+ */
+bool atom::Compare(const atom &ptr)
+{
+  if (nr < ptr.nr)
     return true;
-  } else
+  else
     return false;
 };
 
-/** Output of atom::nr along with all bond partners.
- * \param *AdjacencyFile output stream
- */
-void atom::OutputAdjacency(ofstream *AdjacencyFile) const
-{
-  *AdjacencyFile << nr << "\t";
-  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
-    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
-  *AdjacencyFile << endl;
-};
-
-/** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
- * \param *file output stream
- */
-void atom::OutputOrder(ofstream *file)
-{
-  *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
-  //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
+/** Returns squared distance to a given vector.
+ * \param origin vector to calculate distance to
+ * \return distance squared
+ */
+double atom::DistanceSquaredToVector(Vector &origin)
+{
+  return origin.DistanceSquared(&x);
+};
+
+/** Returns distance to a given vector.
+ * \param origin vector to calculate distance to
+ * \return distance
+ */
+double atom::DistanceToVector(Vector &origin)
+{
+  return origin.Distance(&x);
 };
 
@@ -313,8 +354,75 @@
 };
 
+
+bool operator < (atom &a, atom &b)
+{
+  return a.Compare(b);
+};
+
+/** Output graph info of this atom.
+ * \param *out output stream
+ */
+void GraphNode::OutputGraphInfo(ofstream *out) const
+{
+  *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
+  OutputComponentNumber(out);
+  *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
+};
+
+/** Output a list of flags, stating whether the bond was visited or not.
+ * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
+ * \param *out output stream for debugging
+ */
+void GraphNode::OutputComponentNumber(ofstream *out) const
+{
+  if (ComponentNr != NULL) {
+    for (int i=0; ComponentNr[i] != -1; i++)
+      *out << ComponentNr[i] << " ";
+  }
+};
+
+/** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
+ * \param *file output stream
+ */
+void BondedParticle::OutputOrder(ofstream *file)
+{
+  *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
+  //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
+};
+
+/** Prints all bonds of this atom with total degree.
+ * \param *out stream to output to
+ * \return true - \a *out present, false - \a *out is NULL
+ */
+bool BondedParticle::OutputBondOfAtom(ofstream *out) const
+{
+  if (out != NULL) {
+    *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
+    int TotalDegree = 0;
+    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
+      *out << **Runner << "\t";
+      TotalDegree += (*Runner)->BondDegree;
+    }
+    *out << " -- TotalDegree: " << TotalDegree << endl;
+    return true;
+  } else
+    return false;
+};
+
+/** Output of atom::nr along with all bond partners.
+ * \param *AdjacencyFile output stream
+ */
+void BondedParticle::OutputAdjacency(ofstream *AdjacencyFile) const
+{
+  *AdjacencyFile << nr << "\t";
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
+  *AdjacencyFile << endl;
+};
+
 /** Puts a given bond into atom::ListOfBonds.
  * \param *Binder bond to insert
  */
-bool atom::RegisterBond(bond *Binder)
+bool BondedParticle::RegisterBond(bond *Binder)
 {
   bool status = false;
@@ -335,5 +443,5 @@
  * \param *Binder bond to remove
  */
-bool atom::UnregisterBond(bond *Binder)
+bool BondedParticle::UnregisterBond(bond *Binder)
 {
   bool status = false;
@@ -354,5 +462,5 @@
  * \note Does not do any memory de-allocation.
  */
-void atom::UnregisterAllBond()
+void BondedParticle::UnregisterAllBond()
 {
   ListOfBonds.clear();
@@ -362,5 +470,5 @@
  * \param *out output stream for debugging
  */
-int atom::CorrectBondDegree(ofstream *out)
+int BondedParticle::CorrectBondDegree(ofstream *out)
 {
   int NoBonds = 0;
@@ -395,26 +503,40 @@
 };
 
-ostream & operator << (ostream &ost, const atom &a)
-{
-  ost << "[" << a.Name << "|" << &a << "]";
-  return ost;
-};
-
-ostream & atom::operator << (ostream &ost)
-{
-  ost << "[" << Name << "|" << this << "]";
-  return ost;
-};
-
-/** Compares the indices of \a this atom with a given \a ptr.
- * \param ptr atom to compare index against
- * \return true - this one's is smaller, false - not
- */
-bool atom::Compare(const atom &ptr)
-{
-  if (nr < ptr.nr)
-    return true;
-  else
-    return false;
+/** Adds kinetic energy of this atom to given temperature value.
+ * \param *temperature add on this value
+ * \param step given step of trajectory to add
+ */
+void TrajectoryParticle::AddKineticToTemperature(double *temperature, int step) const
+{
+  for (int i=NDIM;i--;)
+    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
+};
+
+/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
+ * \param startstep trajectory begins at
+ * \param endstep trajectory ends at
+ * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
+ * \param *Force Force matrix to store result in
+ */
+void TrajectoryParticle::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
+{
+  double constant = 10.;
+  TrajectoryParticle *Sprinter = PermutationMap[nr];
+  // set forces
+  for (int i=NDIM;i++;)
+    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
+};
+
+/** Correct velocity against the summed \a CoGVelocity for \a step.
+ * \param *ActualTemp sum up actual temperature meanwhile
+ * \param Step MD step in atom::Tracjetory
+ * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
+ */
+void TrajectoryParticle::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
+{
+  for(int d=0;d<NDIM;d++) {
+    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
+    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
+  }
 };
 
@@ -423,5 +545,5 @@
  * \param MaxSteps
  */
-void atom::ResizeTrajectory(int MaxSteps)
+void TrajectoryParticle::ResizeTrajectory(int MaxSteps)
 {
   if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
@@ -437,5 +559,5 @@
  * \param src index of source step
  */
-void atom::CopyStepOnStep(int dest, int src)
+void TrajectoryParticle::CopyStepOnStep(int dest, int src)
 {
   if (dest == src)  // self assignment check
@@ -455,5 +577,5 @@
  * \param *Force matrix with forces
  */
-void atom::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
+void TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
 {
   //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
@@ -484,70 +606,9 @@
  * \param *TotalVelocity pointer to tota velocity sum
  */
-void atom::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
+void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
 {
   *TotalMass += type->mass;  // sum up total mass
   for(int d=0;d<NDIM;d++) {
     TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
-  }
-};
-
-/** Returns squared distance to a given vector.
- * \param origin vector to calculate distance to
- * \return distance squared
- */
-double atom::DistanceSquaredToVector(Vector &origin)
-{
-  return origin.DistanceSquared(&x);
-};
-
-/** Adds kinetic energy of this atom to given temperature value.
- * \param *temperature add on this value
- * \param step given step of trajectory to add
- */
-void atom::AddKineticToTemperature(double *temperature, int step) const
-{
-  for (int i=NDIM;i--;)
-    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
-};
-
-/** Returns distance to a given vector.
- * \param origin vector to calculate distance to
- * \return distance
- */
-double atom::DistanceToVector(Vector &origin)
-{
-  return origin.Distance(&x);
-};
-
-bool operator < (atom &a, atom &b)
-{
-  return a.Compare(b);
-};
-
-/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
- * \param startstep trajectory begins at
- * \param endstep trajectory ends at
- * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
- * \param *Force Force matrix to store result in
- */
-void atom::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
-{
-  double constant = 10.;
-  atom *Sprinter = PermutationMap[nr];
-  // set forces
-  for (int i=NDIM;i++;)
-    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
-};
-
-/** Correct velocity against the summed \a CoGVelocity for \a step.
- * \param *ActualTemp sum up actual temperature meanwhile
- * \param Step MD step in atom::Tracjetory
- * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
- */
-void atom::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
-{
-  for(int d=0;d<NDIM;d++) {
-    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
-    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
   }
 };
@@ -558,5 +619,5 @@
  * \param *ekin sum of kinetic energy
  */
-void atom::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
+void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
 {
   double *U = Trajectory.U.at(Step).x;
@@ -573,5 +634,5 @@
  * \param *E
  */
-void atom::Thermostat_Gaussian_init(int Step, double *G, double *E)
+void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
 {
   double *U = Trajectory.U.at(Step).x;
@@ -590,5 +651,5 @@
  * \param *configuration configuration class with TempFrequency and TargetTemp
  */
-void atom::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
+void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
 {
   double *U = Trajectory.U.at(Step).x;
@@ -606,5 +667,5 @@
  * \param *configuration configuration class with TempFrequency and TargetTemp
  */
-void atom::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
+void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
 {
   double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
@@ -631,5 +692,5 @@
  * \param *configuration configuration class with TempFrequency and Deltat
  */
-void atom::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
+void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
 {
   double *U = Trajectory.U.at(Step).x;
@@ -646,5 +707,5 @@
  * \param *delta_alpha additional sum of kinetic energy on return
  */
-void atom::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
+void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
 {
   double *U = Trajectory.U.at(Step).x;
@@ -661,5 +722,5 @@
  * \param *configuration configuration class with TempFrequency and Deltat
  */
-void atom::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
+void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
 {
   double *U = Trajectory.U.at(Step).x;
Index: molecuilder/src/atom.hpp
===================================================================
--- molecuilder/src/atom.hpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/atom.hpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -28,4 +28,5 @@
 /****************************************** forward declarations *****************************/
 
+class atom;
 class bond;
 class config;
@@ -38,36 +39,135 @@
 /********************************************** declarations *******************************/
 
+class AtomInfo {
+public:
+  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
+
+  AtomInfo();
+  ~AtomInfo();
+
+private:
+};
+
+class TrajectoryParticleInfo {
+public:
+  struct
+  {
+    vector<Vector> R;  //!< position vector
+    vector<Vector> U;  //!< velocity vector
+    vector<Vector> F;  //!< last force vector
+  } Trajectory;
+  int FixedIon;   //!< config variable that states whether forces act on the ion or not
+
+  TrajectoryParticleInfo();
+  ~TrajectoryParticleInfo();
+
+private:
+};
+
+class TrajectoryParticle :  public TrajectoryParticleInfo, public virtual AtomInfo, public virtual ParticleInfo {
+public:
+
+  TrajectoryParticle();
+  virtual ~TrajectoryParticle();
+
+  // constraint potential and dynamics stuff
+  void AddKineticToTemperature(double *temperature, int step) const;
+  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
+  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
+
+  // trajectory stuff
+  void ResizeTrajectory(int MaxSteps);
+  void CopyStepOnStep(int dest, int src);
+  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
+  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
+
+  // thermostats
+  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
+  void Thermostat_Gaussian_init(int Step, double *G, double *E);
+  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
+  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
+  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
+  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
+  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
+
+private:
+
+};
+
+class GraphNode;
+
+class GraphNodeInfo {
+public:
+
+  int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
+  int *ComponentNr;//!< belongs to this non-separable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
+  int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect non-separable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
+  bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
+  bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
+  atom *Ancestor; //!< "Father" in Depth-First-Search
+
+  GraphNodeInfo();
+  ~GraphNodeInfo();
+private:
+
+};
+
+
+class GraphNode : public GraphNodeInfo, public virtual ParticleInfo {
+public:
+
+  GraphNode();
+  virtual ~GraphNode();
+
+  void OutputGraphInfo(ofstream *out) const;
+  void OutputComponentNumber(ofstream *out) const;
+
+private:
+
+};
+
+class BondedParticleInfo {
+public:
+  unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
+  bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
+  BondList ListOfBonds; //!< list of all bonds
+
+  BondedParticleInfo();
+  ~BondedParticleInfo();
+
+private:
+
+};
+
+class BondedParticle : public BondedParticleInfo, public virtual ParticleInfo, public virtual AtomInfo {
+public:
+  BondedParticle();
+  virtual ~BondedParticle();
+
+  bool RegisterBond(bond *Binder);
+  bool UnregisterBond(bond *Binder);
+  void UnregisterAllBond();
+  int CountBonds() const;
+  int CorrectBondDegree(ofstream *out);
+  bool OutputBondOfAtom(ofstream *out) const;
+  void OutputAdjacency(ofstream *AdjacencyFile) const;
+  void OutputOrder(ofstream *file);
+
+private:
+
+};
+
 /** Single atom.
  * Class incorporates position, type
  */
-class atom : public TesselPoint {
+class atom : public TesselPoint, public TrajectoryParticle, public GraphNode, public BondedParticle, public virtual ParticleInfo, public virtual AtomInfo {
   public:
-    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
-    BondList ListOfBonds; //!< list of all bonds
-    element *type;  //!< pointing to element
     atom *previous; //!< previous atom in molecule list
     atom *next;     //!< next atom in molecule list
     atom *father;   //!< In many-body bond order fragmentations points to originating atom
-    atom *Ancestor; //!< "Father" in Depth-First-Search
-    //char *Name;      //!< unique name used during many-body bond-order fragmentation, comes from TesselPoint
-    int FixedIon;   //!< config variable that states whether forces act on the ion or not
     int *sort;      //!< sort criteria
-    //int nr;         //!< continuous, unique number, comes from TesselPoint
-    int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
-    int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
-    int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
-    bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
-    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
-    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
-    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
 
   atom();
@@ -80,7 +180,7 @@
   bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
   bool OutputTrajectoryXYZ(ofstream *out, int step) const;
-  bool OutputBondOfAtom(ofstream *out) const;
-  void OutputAdjacency(ofstream *AdjacencyFile) const;
+  void OutputMPQCLine(ofstream *out, Vector *center, int *AtomNo) const;
 
+  void InitComponentNr();
 
   void EqualsFather ( atom *ptr, atom **res );
@@ -89,48 +189,10 @@
   bool Compare(const atom &ptr);
 
-  // trajectory stuff
-  void ResizeTrajectory(int MaxSteps);
-  void CopyStepOnStep(int dest, int src);
-  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
-  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
-
   double DistanceToVector(Vector &origin);
   double DistanceSquaredToVector(Vector &origin);
-
   bool IsInParallelepiped(Vector offset, double *parallelepiped);
-
-  // bond order stuff
-  void OutputOrder(ofstream *file);
-  void OutputGraphInfo(ofstream *out) const;
-  void OutputComponentNumber(ofstream *out) const;
-  void InitComponentNr();
-  int CountBonds() const;
-  int CorrectBondDegree(ofstream *out);
-
-  bool RegisterBond(bond *Binder);
-  bool UnregisterBond(bond *Binder);
-  void UnregisterAllBond();
-
-  // constraint potential and dynamics stuff
-  void AddKineticToTemperature(double *temperature, int step) const;
-  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
-  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
-
-  // thermostats
-  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
-  void Thermostat_Gaussian_init(int Step, double *G, double *E);
-  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
-  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
-  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
-  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
-  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
-
-
-  ostream & operator << (ostream &ost);
 
   private:
 };
 
-ostream & operator << (ostream &ost, const atom &a);
-
 #endif /* ATOM_HPP_ */
Index: molecuilder/src/bond.cpp
===================================================================
--- molecuilder/src/bond.cpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/bond.cpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -95,5 +95,5 @@
  * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
  */
-atom * bond::GetOtherAtom(const atom *Atom) const
+atom * bond::GetOtherAtom(const ParticleInfo *Atom) const
 {
   if(leftatom == Atom) 
@@ -118,5 +118,5 @@
  * \return true if it is either bond::leftatom or bond::rightatom, false otherwise
  */
-bool bond::Contains(const atom *ptr)
+bool bond::Contains(const ParticleInfo *ptr)
 {
   return ((leftatom == ptr) || (rightatom == ptr));
Index: molecuilder/src/bond.hpp
===================================================================
--- molecuilder/src/bond.hpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/bond.hpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -21,4 +21,5 @@
 
 class atom;
+class ParticleInfo;
 
 /********************************************** declarations *******************************/
@@ -43,10 +44,10 @@
     enum EdgeType Type;//!< whether this is a tree or back edge
 
-  atom * GetOtherAtom(const atom *Atom) const;
+  atom * GetOtherAtom(const ParticleInfo *Atom) const;
 
   bool MarkUsed(enum Shading color);
   enum Shading IsUsed();
   void ResetUsed();
-  bool Contains(const atom *ptr);
+  bool Contains(const ParticleInfo *ptr);
   bool Contains(const int nr);
 
Index: molecuilder/src/config.cpp
===================================================================
--- molecuilder/src/config.cpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/config.cpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -9,4 +9,5 @@
 #include "element.hpp"
 #include "helpers.hpp"
+#include "molecule.hpp"
 #include "memoryallocator.hpp"
 #include "molecule.hpp"
@@ -1444,8 +1445,5 @@
 bool config::SaveMPQC(const char *filename, molecule *mol) const
 {
-  int ElementNo = 0;
   int AtomNo;
-  atom *Walker = NULL;
-  element *runner = NULL;
   Vector *center = NULL;
   ofstream *output = NULL;
@@ -1475,20 +1473,6 @@
   center = mol->DetermineCenterOfAll(output);
   // output of atoms
-  runner = mol->elemente->start;
-  while (runner->next != mol->elemente->end) { // go through every element
-    runner = runner->next;
-    if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
-      ElementNo++;
-      AtomNo = 0;
-      Walker = mol->start;
-      while (Walker->next != mol->end) { // go through every atom of this element
-        Walker = Walker->next;
-        if (Walker->type == runner) { // if this atom fits to element
-          AtomNo++;
-          *output << "\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
-        }
-      }
-    }
-  }
+  AtomNo = 0;
+  mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, center, &AtomNo );
   delete(center);
   *output << "\t}" << endl;
@@ -1524,20 +1508,6 @@
   center = mol->DetermineCenterOfAll(output);
   // output of atoms
-  runner = mol->elemente->start;
-  while (runner->next != mol->elemente->end) { // go through every element
-    runner = runner->next;
-    if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
-      ElementNo++;
-      AtomNo = 0;
-      Walker = mol->start;
-      while (Walker->next != mol->end) { // go through every atom of this element
-        Walker = Walker->next;
-        if (Walker->type == runner) { // if this atom fits to element
-          AtomNo++;
-          *output << "\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
-        }
-      }
-    }
-  }
+  AtomNo = 0;
+  mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, center, &AtomNo );
   delete(center);
   *output << "\t}" << endl;
Index: molecuilder/src/molecule.cpp
===================================================================
--- molecuilder/src/molecule.cpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/molecule.cpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -1078,5 +1078,5 @@
   for (int step=startstep;step < endstep; step++) { // loop over all time steps
     temperature = 0.;
-    ActOnAllAtoms( &atom::AddKineticToTemperature, &temperature, step);
+    ActOnAllAtoms( &TrajectoryParticle::AddKineticToTemperature, &temperature, step);
     *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
   }
@@ -1084,5 +1084,5 @@
 };
 
-void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int TesselPoint::*index)
+void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index)
 {
   atom *Walker = start;
Index: molecuilder/src/molecule.hpp
===================================================================
--- molecuilder/src/molecule.hpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/molecule.hpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -146,24 +146,24 @@
 
   // templates for allowing global manipulation of all atoms
-  template <typename res> void ActOnAllAtoms( res (atom::*f)() );
-  template <typename res> void ActOnAllAtoms( res (atom::*f)() const );
-  template <typename res> void ActOnAllAtoms( res (atom::*f)() ) const;
-  template <typename res> void ActOnAllAtoms( res (atom::*f)() const) const;
-  template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T), T t );
-  template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T) const, T t );
-  template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T), T t ) const;
-  template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T) const, T t ) const;
-  template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U), T t, U u );
-  template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u );
-  template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U), T t, U u ) const;
-  template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u ) const;
-  template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v);
-  template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V) const, T t, U u, V v);
-  template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v) const;
-  template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V) const, T t, U u, V v) const;
-  template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W), T t, U u, V v, W w);
-  template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W) const, T t, U u, V v, W w);
-  template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W), T t, U u, V v, W w) const;
-  template <typename res, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (atom::*f)(T, U, V, W) const, T t, U u, V v, W w) const;
+  template <typename res, typename typ> void ActOnAllAtoms( res (typ::*f)() );
+  template <typename res, typename typ> void ActOnAllAtoms( res (typ::*f)() const );
+  template <typename res, typename typ> void ActOnAllAtoms( res (typ::*f)() ) const;
+  template <typename res, typename typ> void ActOnAllAtoms( res (typ::*f)() const) const;
+  template <typename res, typename typ, typename T> void ActOnAllAtoms( res (typ::*f)(T), T t );
+  template <typename res, typename typ, typename T> void ActOnAllAtoms( res (typ::*f)(T) const, T t );
+  template <typename res, typename typ, typename T> void ActOnAllAtoms( res (typ::*f)(T), T t ) const;
+  template <typename res, typename typ, typename T> void ActOnAllAtoms( res (typ::*f)(T) const, T t ) const;
+  template <typename res, typename typ, typename T, typename U> void ActOnAllAtoms( res (typ::*f)(T, U), T t, U u );
+  template <typename res, typename typ, typename T, typename U> void ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u );
+  template <typename res, typename typ, typename T, typename U> void ActOnAllAtoms( res (typ::*f)(T, U), T t, U u ) const;
+  template <typename res, typename typ, typename T, typename U> void ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u ) const;
+  template <typename res, typename typ, typename T, typename U, typename V> void ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v);
+  template <typename res, typename typ, typename T, typename U, typename V> void ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v);
+  template <typename res, typename typ, typename T, typename U, typename V> void ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v) const;
+  template <typename res, typename typ, typename T, typename U, typename V> void ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v) const;
+  template <typename res, typename typ, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w);
+  template <typename res, typename typ, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w);
+  template <typename res, typename typ, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w) const;
+  template <typename res, typename typ, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w) const;
 
   // templates for allowing conditional global copying of molecule with each atom as single argument
@@ -174,27 +174,26 @@
 
   // templates for allowing global manipulation of an array with one entry per atom
-  void SetIndexedArrayForEachAtomTo ( atom **array, int TesselPoint::* index);
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T *));
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T *), T t);
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T *), T *t);
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *));
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T t);
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T *t);
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, T (atom::*Setor)(Vector &), Vector atom::*value);
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, T (atom::*Setor)(Vector &), Vector &vect );
+  void SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::* index);
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *));
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *), T t);
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *), T *t);
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *));
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *), T t);
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *), T *t);
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(Vector &), Vector atom::*value);
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(Vector &), Vector &vect );
 
   // templates for allowing global manipulation of each atom by entries in an array
-  template <typename T> void SetAtomValueToIndexedArray ( T *array, int TesselPoint::*index, T atom::*value );
-  template <typename T> void SetAtomValueToIndexedArray ( T *array, int element::*index, T atom::*value );
-  template <typename T> void SetAtomValueToValue ( T value, T atom::*ptr );
-
-  template <typename res> res SumPerAtom(res (atom::*f)() );
-  template <typename res> res SumPerAtom(res (atom::*f)() const );
-  template <typename res> res SumPerAtom(res (atom::*f)() ) const;
-  template <typename res> res SumPerAtom(res (atom::*f)() const ) const;
-  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) , T t );
-  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) const, T t );
-  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) , T t ) const;
-  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) const, T t ) const;
+  template <typename T, typename typ, typename typ2> void SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value );
+  template <typename T, typename typ> void SetAtomValueToValue ( T value, T typ::*ptr );
+
+  template <typename res, typename typ> res SumPerAtom(res (typ::*f)() );
+  template <typename res, typename typ> res SumPerAtom(res (typ::*f)() const );
+  template <typename res, typename typ> res SumPerAtom(res (typ::*f)() ) const;
+  template <typename res, typename typ> res SumPerAtom(res (typ::*f)() const ) const;
+  template <typename res, typename typ, typename T> res SumPerAtom(res (typ::*f)(T) , T t );
+  template <typename res, typename typ, typename T> res SumPerAtom(res (typ::*f)(T) const, T t );
+  template <typename res, typename typ, typename T> res SumPerAtom(res (typ::*f)(T) , T t ) const;
+  template <typename res, typename typ, typename T> res SumPerAtom(res (typ::*f)(T) const, T t ) const;
 
   /// remove atoms from molecule.
Index: molecuilder/src/molecule_template.hpp
===================================================================
--- molecuilder/src/molecule_template.hpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/molecule_template.hpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -156,5 +156,5 @@
 
 // zero arguments
-template <typename res> res molecule::SumPerAtom(res (atom::*f)() )
+template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() )
 {
   res result = 0;
@@ -166,5 +166,5 @@
   return result;
 };
-template <typename res> res molecule::SumPerAtom(res (atom::*f)() const )
+template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() const )
 {
   res result = 0;
@@ -176,5 +176,5 @@
   return result;
 };
-template <typename res> res molecule::SumPerAtom(res (atom::*f)() ) const
+template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() ) const
 {
   res result = 0;
@@ -186,5 +186,5 @@
   return result;
 };
-template <typename res> res molecule::SumPerAtom(res (atom::*f)() const ) const
+template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() const ) const
 {
   res result = 0;
@@ -197,5 +197,5 @@
 };
 // one argument
-template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T), T t )
+template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T), T t )
 {
   res result = 0;
@@ -207,5 +207,5 @@
   return result;
 };
-template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T) const, T t )
+template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T) const, T t )
 {
   res result = 0;
@@ -217,5 +217,5 @@
   return result;
 };
-template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T), T t ) const
+template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T), T t ) const
 {
   res result = 0;
@@ -227,5 +227,5 @@
   return result;
 };
-template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T) const, T t ) const
+template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T) const, T t ) const
 {
   res result = 0;
@@ -357,5 +357,5 @@
 
 // zero arguments
-template <typename res> void molecule::ActOnAllAtoms( res (atom::*f)())
+template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)())
 {
   atom *Walker = start;
@@ -365,5 +365,5 @@
   }
 };
-template <typename res> void molecule::ActOnAllAtoms( res (atom::*f)() const)
+template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const)
 {
   atom *Walker = start;
@@ -373,5 +373,5 @@
   }
 };
-template <typename res> void molecule::ActOnAllAtoms( res (atom::*f)()) const
+template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)()) const
 {
   atom *Walker = start;
@@ -381,5 +381,5 @@
   }
 };
-template <typename res> void molecule::ActOnAllAtoms( res (atom::*f)() const) const
+template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const) const
 {
   atom *Walker = start;
@@ -390,5 +390,5 @@
 };
 // one argument
-template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T), T t )
+template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t )
 {
   atom *Walker = start;
@@ -398,5 +398,5 @@
   }
 };
-template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T) const, T t )
+template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t )
 {
   atom *Walker = start;
@@ -406,5 +406,5 @@
   }
 };
-template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T), T t ) const
+template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t ) const
 {
   atom *Walker = start;
@@ -414,5 +414,5 @@
   }
 };
-template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T) const, T t ) const
+template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t ) const
 {
   atom *Walker = start;
@@ -423,5 +423,5 @@
 };
 // two argument
-template <typename res, typename T, typename U> void molecule::ActOnAllAtoms( res (atom::*f)(T, U), T t, U u )
+template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U), T t, U u )
 {
   atom *Walker = start;
@@ -431,5 +431,5 @@
   }
 };
-template <typename res, typename T, typename U> void molecule::ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u )
+template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u )
 {
   atom *Walker = start;
@@ -439,5 +439,5 @@
   }
 };
-template <typename res, typename T, typename U> void molecule::ActOnAllAtoms( res (atom::*f)(T, U), T t, U u ) const
+template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U), T t, U u ) const
 {
   atom *Walker = start;
@@ -447,5 +447,5 @@
   }
 };
-template <typename res, typename T, typename U> void molecule::ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u ) const
+template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u ) const
 {
   atom *Walker = start;
@@ -456,5 +456,5 @@
 };
 // three argument
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v)
+template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v)
 {
   atom *Walker = start;
@@ -464,5 +464,5 @@
   }
 };
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V) const, T t, U u, V v)
+template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v)
 {
   atom *Walker = start;
@@ -472,5 +472,5 @@
   }
 };
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v) const
+template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v) const
 {
   atom *Walker = start;
@@ -480,5 +480,5 @@
   }
 };
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V) const, T t, U u, V v) const
+template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v) const
 {
   atom *Walker = start;
@@ -489,5 +489,5 @@
 };
 // four arguments
-template <typename res, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V, W), T t, U u, V v, W w)
+template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w)
 {
   atom *Walker = start;
@@ -497,5 +497,5 @@
   }
 };
-template <typename res, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V, W) const, T t, U u, V v, W w)
+template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w)
 {
   atom *Walker = start;
@@ -505,5 +505,5 @@
   }
 };
-template <typename res, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V, W), T t, U u, V v, W w) const
+template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w) const
 {
   atom *Walker = start;
@@ -513,5 +513,5 @@
   }
 };
-template <typename res, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V, W) const, T t, U u, V v, W w) const
+template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w) const
 {
   atom *Walker = start;
@@ -525,5 +525,5 @@
 
 // for atom ints
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, void (*Setor)(T *, T *) )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *) )
 {
   atom *Walker = start;
@@ -534,5 +534,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, void (*Setor)(T *, T *), T value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T value )
 {
   atom *Walker = start;
@@ -542,5 +542,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, void (*Setor)(T *, T *), T *value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T *value )
 {
   atom *Walker = start;
@@ -550,5 +550,5 @@
   }
 };
-// for element ints inside atom class
+// for element ints
 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *) )
 {
@@ -576,5 +576,6 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, T (atom::*Setor)(Vector &), Vector atom::*value )
+
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(Vector &), Vector atom::*value )
 {
   atom *Walker = start;
@@ -584,5 +585,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, T (atom::*Setor)(Vector &), Vector &vect )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(Vector &), Vector &vect )
 {
   atom *Walker = start;
@@ -593,5 +594,5 @@
 };
 
-template <typename T> void molecule::SetAtomValueToIndexedArray ( T *array, int TesselPoint::*index, T atom::*value )
+template <typename T, typename typ, typename typ2> void molecule::SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value )
 {
   atom *Walker = start;
@@ -603,14 +604,5 @@
 };
 
-template <typename T> void molecule::SetAtomValueToIndexedArray ( T *array, int element::*index, T atom::*value )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    Walker->*value = array[(Walker->type->*index)];
-    //cout << Verbose(2) << *Walker << " gets " << (Walker->*value) << endl;
-  }
-};
-template <typename T> void molecule::SetAtomValueToValue ( T value, T atom::*ptr )
+template <typename T, typename typ> void molecule::SetAtomValueToValue ( T value, T typ::*ptr )
 {
   atom *Walker = start;
Index: molecuilder/src/tesselation.cpp
===================================================================
--- molecuilder/src/tesselation.cpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/tesselation.cpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -16,4 +16,17 @@
 
 class molecule;
+
+// =========================================== ParticleInfo ====================================
+
+/** Constructor of ParticleInfo.
+ */
+ParticleInfo::ParticleInfo() : nr(-1), Name(NULL) {};
+
+/** Destructor of ParticleInfo.
+ */
+ParticleInfo::~ParticleInfo()
+{
+  Free(&Name);
+};
 
 // ======================================== Points on Boundary =================================
Index: molecuilder/src/tesselation.hpp
===================================================================
--- molecuilder/src/tesselation.hpp	(revision 99bed31cf29f533f16d3526fc03ba2d88887bf67)
+++ molecuilder/src/tesselation.hpp	(revision 6b937bdc52fc2fbea7d288c718973ad26f2fda1b)
@@ -25,4 +25,5 @@
 
 #include "vector.hpp"
+#include "helpers.hpp"
 
 /****************************************** forward declarations *****************************/
@@ -144,7 +145,22 @@
 // =========================================================== class TESSELPOINT ===========================================
 
+class ParticleInfo {
+public:
+  int nr;       // index to easierly identify
+  char *Name;   // some name to reference to on output
+
+  ParticleInfo();
+  ~ParticleInfo();
+
+  ostream & operator << (ostream &ost);
+
+private:
+};
+
+ostream & operator << (ostream &ost, const ParticleInfo &a);
+
 /** Is a single point of the set of Vectors, also a super-class to be inherited and and its functions to be implemented.
  */
-class TesselPoint {
+class TesselPoint : virtual public ParticleInfo {
 public:
   TesselPoint();
@@ -152,6 +168,4 @@
 
   Vector *node;   // pointer to position of the dot in space
-  int nr;       // index to easierly identify
-  char *Name;   // some name to reference to on output
 
   virtual ostream & operator << (ostream &ost);
