Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision f66195d268aa8acd6e1a3d8ddc7ba7f220eb784b)
+++ src/atom.cpp	(revision 49f802cfc32887da6d23e51283966873ce964f38)
@@ -222,11 +222,17 @@
 {
   if (out != NULL) {
-    *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: ";
-    int TotalDegree = 0;
-    for (int j=0;j<NumberOfBondsPerAtom[nr];j++) {
-      *out << *ListOfBondsPerAtom[nr][j] << "\t";
-      TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree;
+#ifdef ADDHYDROGEN
+    if (type->Z != 1) {   // regard only non-hydrogen
+#endif
+      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: ";
+      int TotalDegree = 0;
+      for (int j=0;j<NumberOfBondsPerAtom[nr];j++) {
+        *out << *ListOfBondsPerAtom[nr][j] << "\t";
+        TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree;
+      }
+      *out << " -- TotalDegree: " << TotalDegree << endl;
+#ifdef ADDHYDROGEN
     }
-    *out << " -- TotalDegree: " << TotalDegree << endl;
+#endif
     return true;
   } else
@@ -258,4 +264,32 @@
 };
 
+/** 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)
 {
Index: src/atom.hpp
===================================================================
--- src/atom.hpp	(revision f66195d268aa8acd6e1a3d8ddc7ba7f220eb784b)
+++ src/atom.hpp	(revision 49f802cfc32887da6d23e51283966873ce964f38)
@@ -79,4 +79,9 @@
   bool Compare(const atom &ptr);
 
+  double DistanceToVector(Vector &origin);
+  double DistanceSquaredToVector(Vector &origin);
+
+  void AddKineticToTemperature(double *temperature, int step) const;
+
   bool IsInParallelepiped(Vector offset, double *parallelepiped);
 
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision f66195d268aa8acd6e1a3d8ddc7ba7f220eb784b)
+++ src/molecule.cpp	(revision 49f802cfc32887da6d23e51283966873ce964f38)
@@ -779,18 +779,5 @@
 {
   *out << Verbose(2) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-#ifdef ADDHYDROGEN
-    if (Walker->type->Z != 1) {   // regard only non-hydrogen
-#endif
-      *out << Verbose(2) << "Atom " << Walker->Name << " has Bonds: "<<endl;
-      for(int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {
-        *out << Verbose(3) << *(ListOfBondsPerAtom)[Walker->nr][j] << endl;
-      }
-#ifdef ADDHYDROGEN
-    }
-#endif
-  }
+  ActOnAllAtoms (&atom::OutputBondOfAtom, out, NumberOfBondsPerAtom, ListOfBondsPerAtom);
   *out << endl;
 };
@@ -983,5 +970,4 @@
   size_t *PermMap = NULL, *OtherPermMap = NULL;
   int *PermutationMap = NULL;
-  atom *Walker = NULL;
   bool result = true; // status of comparison
 
@@ -1042,14 +1028,6 @@
     Distances = Malloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
     OtherDistances = Malloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      Distances[Walker->nr] = CenterOfGravity.DistanceSquared(&Walker->x);
-    }
-    Walker = OtherMolecule->start;
-    while (Walker->next != OtherMolecule->end) {
-      Walker = Walker->next;
-      OtherDistances[Walker->nr] = OtherCenterOfGravity.DistanceSquared(&Walker->x);
-    }
+    SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, CenterOfGravity);
+    SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, CenterOfGravity);
 
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
@@ -1150,5 +1128,4 @@
 {
   double temperature;
-  atom *Walker = NULL;
   // test stream
   if (output == NULL)
@@ -1158,10 +1135,5 @@
   for (int step=startstep;step < endstep; step++) { // loop over all time steps
     temperature = 0.;
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      for (int i=NDIM;i--;)
-        temperature += Walker->type->mass * Walker->Trajectory.U.at(step).x[i]* Walker->Trajectory.U.at(step).x[i];
-    }
+    ActOnAllAtoms( &atom::AddKineticToTemperature, &temperature, step);
     *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
   }
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision f66195d268aa8acd6e1a3d8ddc7ba7f220eb784b)
+++ src/molecule.hpp	(revision 49f802cfc32887da6d23e51283966873ce964f38)
@@ -108,13 +108,17 @@
   template <typename res> void ActOnAllVectors( res (Vector::*f)() const);
   template <typename res> void ActOnAllVectors( res (Vector::*f)() ) const;
+  template <typename res> void ActOnAllVectors( res (Vector::*f)() const) const;
   template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t );
   template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t );
   template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t ) const;
+  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t ) const;
   template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u );
   template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u );
   template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const;
+  template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const;
   template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v);
   template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v);
   template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const;
+  template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const;
 
   // templates for allowing global manipulation of molecule with each atom as single argument
@@ -122,4 +126,5 @@
   template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) const);
   template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) ) const;
+  template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) const) const;
 
   // templates for allowing global copying of molecule with each atom as single argument
@@ -127,4 +132,5 @@
   template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const , molecule *copy);
   template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const;
+  template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const;
 
   // templates for allowing global manipulation of all atoms
@@ -132,16 +138,21 @@
   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;
 
   // templates for allowing conditional global copying of molecule with each atom as single argument
@@ -152,8 +163,10 @@
 
   // templates for allowing global manipulation of an array with one entry per atom
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int atom::* index, void (*Setor)(T *, T));
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int atom::* index, void (*Setor)(T *, T), T);
+  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 element::*index, void (*Setor)(T *, T));
-  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(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 );
 
   /// remove atoms from molecule.
Index: src/molecule_template.hpp
===================================================================
--- src/molecule_template.hpp	(revision f66195d268aa8acd6e1a3d8ddc7ba7f220eb784b)
+++ src/molecule_template.hpp	(revision 49f802cfc32887da6d23e51283966873ce964f38)
@@ -45,4 +45,12 @@
   }
 };
+template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() const ) const
+    {
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    ((Walker->node)->*f)();
+  }
+};
 // one argument
 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t )
@@ -70,4 +78,12 @@
   }
 };
+template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    ((Walker->node)->*f)(t);
+  }
+};
 // two arguments
 template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u )
@@ -95,4 +111,12 @@
   }
 };
+template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    ((Walker->node)->*f)(t, u);
+  }
+};
 // three arguments
 template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v)
@@ -120,4 +144,12 @@
   }
 };
+template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    ((Walker->node)->*f)(t, u, v);
+  }
+};
 
 // ================== Acting with each Atoms on same molecule ========================== //
@@ -148,4 +180,12 @@
   }
 };
+template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *) const) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (*f)(Walker);
+  }
+};
 
 // ================== Acting with each Atoms on copy molecule ========================== //
@@ -176,4 +216,12 @@
   }
 };
+template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (copy->*f)(Walker);
+  }
+};
 
 // ================== Acting with each Atoms on copy molecule if true ========================== //
@@ -247,4 +295,12 @@
   }
 };
+template <typename res> void molecule::ActOnAllAtoms( res (atom::*f)() const) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)();
+  }
+};
 // one argument
 template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T), T t )
@@ -264,4 +320,12 @@
   }
 };
+template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T), T t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)(t);
+  }
+};
 template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T) const, T t ) const
 {
@@ -297,4 +361,12 @@
   }
 };
+template <typename res, typename T, typename U> void molecule::ActOnAllAtoms( res (atom::*f)(T, U) const, T t, U u ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)(t, u);
+  }
+};
 // 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)
@@ -322,4 +394,12 @@
   }
 };
+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
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)(t, u, v);
+  }
+};
 // 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)
@@ -347,9 +427,17 @@
   }
 };
+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
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)(t, u, v, w);
+  }
+};
 
 // ===================== Accesssing arrays indexed by some integer for each atom ======================
 
 // for atom ints
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int atom::*index, void (*Setor)(T *, T) )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, void (*Setor)(T *, T) )
 {
   atom *Walker = start;
@@ -359,5 +447,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int atom::*index, void (*Setor)(T *, T), T value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, void (*Setor)(T *, T), T value )
 {
   atom *Walker = start;
@@ -384,4 +472,20 @@
   }
 };
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, T (atom::*Setor)(Vector &), Vector atom::*value )
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    array[(Walker->*index)] = (Walker->*Setor) (Walker->*value);
+  }
+};
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::*index, T (atom::*Setor)(Vector &), Vector &vect )
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    array[(Walker->*index)] = (Walker->*Setor) (vect);
+  }
+};
 
 #endif /* MOLECULE_TEMPLATE_HPP_ */
