Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -43,5 +43,5 @@
   sort = &nr;
   node = &x;
-}
+};
 
 
@@ -83,5 +83,5 @@
  * \param **res return value (only set if atom::father is equal to \a *ptr)
  */
-void atom::EqualsFather ( atom *ptr, atom **res )
+void atom::EqualsFather ( const atom *ptr, const atom **res ) const
 {
   if ( ptr == father )
@@ -94,5 +94,5 @@
  * \return true - is inside, false - is not
  */
-bool atom::IsInParallelepiped(Vector offset, double *parallelepiped)
+bool atom::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
 {
   return (node->IsInParallelepiped(offset, parallelepiped));
@@ -110,5 +110,5 @@
 };
 
-/** Output of a single atom.
+/** Output of a single atom with given numbering.
  * \param ElementNo cardinal number of the element
  * \param AtomNo cardinal number among these atoms of the same element
@@ -117,5 +117,5 @@
   * \return true - \a *out present, false - \a *out is NULL
  */
-bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const
+bool atom::OutputIndexed(ofstream *out, const int ElementNo, const int AtomNo, const char *comment) const
 {
   if (out != NULL) {
@@ -133,5 +133,13 @@
     return false;
 };
-bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
+
+/** Output of a single atom with numbering from array according to atom::type.
+ * \param *ElementNo cardinal number of the element
+ * \param *AtomNo cardinal number among these atoms of the same element
+ * \param *out stream to output to
+ * \param *comment commentary after '#' sign
+  * \return true - \a *out present, false - \a *out is NULL
+ */
+bool atom::OutputArrayIndexed(ofstream *out, const int *ElementNo, int *AtomNo, const char *comment) const
 {
   AtomNo[type->Z]++;  // increment number
@@ -171,5 +179,5 @@
   * \return true - \a *out present, false - \a *out is NULL
  */
-bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const
+bool atom::OutputTrajectory(ofstream *out, const int *ElementNo, int *AtomNo, const int step) const
 {
   AtomNo[type->Z]++;
@@ -193,5 +201,5 @@
  * \return true - \a *out present, false - \a *out is NULL
  */
-bool atom::OutputTrajectoryXYZ(ofstream *out, int step) const
+bool atom::OutputTrajectoryXYZ(ofstream *out, const int step) const
 {
   if (out != NULL) {
@@ -210,5 +218,5 @@
  * \param *AtomNo pointer to atom counter that is increased by one
  */
-void atom::OutputMPQCLine(ofstream *out, Vector *center, int *AtomNo = NULL) const
+void atom::OutputMPQCLine(ofstream *out, const 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;
@@ -217,21 +225,9 @@
 };
 
-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)
+bool atom::Compare(const atom &ptr) const
 {
   if (nr < ptr.nr)
@@ -245,5 +241,5 @@
  * \return distance squared
  */
-double atom::DistanceSquaredToVector(Vector &origin)
+double atom::DistanceSquaredToVector(const Vector &origin) const
 {
   return origin.DistanceSquared(&x);
@@ -254,5 +250,5 @@
  * \return distance
  */
-double atom::DistanceToVector(Vector &origin)
+double atom::DistanceToVector(const Vector &origin) const
 {
   return origin.Distance(&x);
Index: src/atom.hpp
===================================================================
--- src/atom.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -49,21 +49,21 @@
   virtual ~atom();
 
-  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 OutputIndexed(ofstream *out, const int ElementNo, const int AtomNo, const char *comment = NULL) const;
+  bool OutputArrayIndexed(ofstream *out, const int *ElementNo, int *AtomNo, const char *comment = NULL) const;
   bool OutputXYZLine(ofstream *out) const;
-  bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
-  bool OutputTrajectoryXYZ(ofstream *out, int step) const;
-  void OutputMPQCLine(ofstream *out, Vector *center, int *AtomNo) const;
+  bool OutputTrajectory(ofstream *out, const int *ElementNo, int *AtomNo, const int step) const;
+  bool OutputTrajectoryXYZ(ofstream *out, const int step) const;
+  void OutputMPQCLine(ofstream *out, const Vector *center, int *AtomNo) const;
 
   void InitComponentNr();
 
-  void EqualsFather ( atom *ptr, atom **res );
+  void EqualsFather ( const atom *ptr, const atom **res ) const;
   void CorrectFather();
   atom *GetTrueFather();
-  bool Compare(const atom &ptr);
+  bool Compare(const atom &ptr) const;
 
-  double DistanceToVector(Vector &origin);
-  double DistanceSquaredToVector(Vector &origin);
-  bool IsInParallelepiped(Vector offset, double *parallelepiped);
+  double DistanceToVector(const Vector &origin) const;
+  double DistanceSquaredToVector(const Vector &origin) const;
+  bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
 
   private:
Index: src/atom_bondedparticle.cpp
===================================================================
--- src/atom_bondedparticle.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom_bondedparticle.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -31,5 +31,5 @@
  * \param *file output stream
  */
-void BondedParticle::OutputOrder(ofstream *file)
+void BondedParticle::OutputOrder(ofstream *file) const
 {
   *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
Index: src/atom_bondedparticle.hpp
===================================================================
--- src/atom_bondedparticle.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom_bondedparticle.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -44,5 +44,5 @@
   bool OutputBondOfAtom(ofstream *out) const;
   void OutputAdjacency(ofstream *AdjacencyFile) const;
-  void OutputOrder(ofstream *file);
+  void OutputOrder(ofstream *file) const;
 
 private:
Index: src/atom_particleinfo.cpp
===================================================================
--- src/atom_particleinfo.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom_particleinfo.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -20,2 +20,14 @@
 };
 
+ostream & operator << (ostream &ost, const ParticleInfo &a)
+{
+  ost << "[" << a.Name << "|" << &a << "]";
+  return ost;
+};
+
+ostream & ParticleInfo::operator << (ostream &ost) const
+{
+  ost << "[" << Name << "|" << this << "]";
+  return ost;
+};
+
Index: src/atom_particleinfo.hpp
===================================================================
--- src/atom_particleinfo.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom_particleinfo.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -33,5 +33,5 @@
   ~ParticleInfo();
 
-  ostream & operator << (ostream &ost);
+  ostream & operator << (ostream &ost) const;
 
 private:
Index: src/atom_trajectoryparticle.cpp
===================================================================
--- src/atom_trajectoryparticle.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom_trajectoryparticle.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -42,5 +42,5 @@
  * \param *Force Force matrix to store result in
  */
-void TrajectoryParticle::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
+void TrajectoryParticle::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force) const
 {
   double constant = 10.;
@@ -129,5 +129,5 @@
  * \param *TotalVelocity pointer to tota velocity sum
  */
-void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
+void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity ) const
 {
   *TotalMass += type->mass;  // sum up total mass
Index: src/atom_trajectoryparticle.hpp
===================================================================
--- src/atom_trajectoryparticle.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/atom_trajectoryparticle.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -44,5 +44,5 @@
   // constraint potential and dynamics stuff
   void AddKineticToTemperature(double *temperature, int step) const;
-  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
+  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force) const;
   void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
 
@@ -51,5 +51,5 @@
   void CopyStepOnStep(int dest, int src);
   void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
-  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
+  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity ) const;
 
   // thermostats
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/config.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -1474,5 +1474,5 @@
   // output of atoms
   AtomNo = 0;
-  mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, center, &AtomNo );
+  mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
   delete(center);
   *output << "\t}" << endl;
@@ -1509,5 +1509,5 @@
   // output of atoms
   AtomNo = 0;
-  mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, center, &AtomNo );
+  mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
   delete(center);
   *output << "\t}" << endl;
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/molecule.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -497,6 +497,6 @@
 
     // get the pendant atoms of current bond in the copy molecule
-    copy->ActOnAllAtoms( &atom::EqualsFather, Binder->leftatom, &LeftAtom );
-    copy->ActOnAllAtoms( &atom::EqualsFather, Binder->rightatom, &RightAtom );
+    copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->leftatom, (const atom **)&LeftAtom );
+    copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->rightatom, (const atom **)&RightAtom );
 
     NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
@@ -526,5 +526,5 @@
  * @param three vectors forming the matrix that defines the shape of the parallelpiped
  */
-molecule* molecule::CopyMoleculeFromSubRegion(Vector offset, double *parallelepiped) {
+molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
   molecule *copy = new molecule(elemente);
 
@@ -733,5 +733,5 @@
         ElementNo[i] = current++;
     }
-    ActOnAllAtoms( &atom::Output, out, ElementNo, AtomNo, (const char *) NULL ); // (bool (atom::*)(int *, int *, ofstream *, const char *))
+    ActOnAllAtoms( &atom::OutputArrayIndexed, (ofstream *)out, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
     return true;
   }
@@ -765,5 +765,5 @@
           ElementNo[i] = current++;
       }
-      ActOnAllAtoms( &atom::OutputTrajectory, out, ElementNo, AtomNo, step );
+      ActOnAllAtoms( &atom::OutputTrajectory, out, (const int *)ElementNo, AtomNo, (const int)step );
     }
     return true;
@@ -971,6 +971,6 @@
     Distances = Malloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
     OtherDistances = Malloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
-    SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, CenterOfGravity);
-    SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, CenterOfGravity);
+    SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
+    SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
 
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
@@ -1084,5 +1084,5 @@
 };
 
-void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index)
+void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
 {
   atom *Walker = start;
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/molecule.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -116,82 +116,70 @@
 
   // templates for allowing global manipulation of all vectors
-  template <typename res> void ActOnAllVectors( res (Vector::*f)() );
-  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
-  template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) );
-  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
-  template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy);
-  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
-  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
-  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () );
-  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t );
-  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u );
-  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v );
+  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () ) const;
+  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () const ) const;
+  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () ) const;
+  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () const ) const;
+  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t ) const;
+  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T) const, T t ) const;
+  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T), T t ) const;
+  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T) const, T t ) const;
+  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const;
+  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const;
+  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const;
+  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const;
+  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const;
+  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const;
+  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const;
+  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const;
 
   // templates for allowing global manipulation of an array with one entry per atom
-  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 );
+  void SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::* index) const;
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *)) const;
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *), T t) const;
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *), T *t) const;
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *)) const;
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *), T t) const;
+  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *), T *t) const;
+  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ atom::*value) const;
+  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ atom::*value) const;
+  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ &vect ) const;
+  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ &vect ) const;
 
   // templates for allowing global manipulation of each atom by entries in an array
-  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 T, typename typ, typename typ2> void SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value ) const;
+  template <typename T, typename typ> void SetAtomValueToValue ( T value, T typ::*ptr ) 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;
@@ -276,5 +264,5 @@
 
   molecule *CopyMolecule();
-  molecule* CopyMoleculeFromSubRegion(Vector offset, double *parallelepiped);
+  molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;
 
   /// Fragment molecule by two different approaches:
Index: src/molecule_template.hpp
===================================================================
--- src/molecule_template.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/molecule_template.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -21,20 +21,4 @@
 
 // zero arguments
-template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)();
-  }
-};
-template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() const )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)();
-  }
-};
 template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() ) const
     {
@@ -54,5 +38,5 @@
 };
 // one argument
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t )
+template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t ) const
 {
   atom *Walker = start;
@@ -62,5 +46,5 @@
   }
 };
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t )
+template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t ) const
 {
   atom *Walker = start;
@@ -70,22 +54,6 @@
   }
 };
-template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t ) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t);
-  }
-};
-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 )
+template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const
 {
   atom *Walker = start;
@@ -95,5 +63,5 @@
   }
 };
-template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u )
+template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const
 {
   atom *Walker = start;
@@ -103,22 +71,6 @@
   }
 };
-template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t, u);
-  }
-};
-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)
+template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const
 {
   atom *Walker = start;
@@ -128,5 +80,5 @@
   }
 };
-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)
+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;
@@ -136,25 +88,9 @@
   }
 };
-template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    ((Walker->node)->*f)(t, u, v);
-  }
-};
-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);
-  }
-};
 
 // ========================= Summing over each Atoms =================================== //
 
 // zero arguments
-template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() )
+template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() ) const
 {
   res result = 0;
@@ -166,5 +102,5 @@
   return result;
 };
-template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() const )
+template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() const ) const
 {
   res result = 0;
@@ -176,5 +112,6 @@
   return result;
 };
-template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() ) const
+// one argument
+template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T), T t ) const
 {
   res result = 0;
@@ -182,9 +119,9 @@
   while (Walker->next != end) {
     Walker = Walker->next;
-    result += (Walker->*f)();
+    result += (Walker->*f)(t);
   }
   return result;
 };
-template <typename res, typename typ> res molecule::SumPerAtom(res (typ::*f)() const ) const
+template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T) const, T t ) const
 {
   res result = 0;
@@ -192,203 +129,222 @@
   while (Walker->next != end) {
     Walker = Walker->next;
-    result += (Walker->*f)();
+    result += (Walker->*f)(t);
   }
   return result;
 };
+
+
+// ================== Acting with each Atoms on same molecule ========================== //
+
+// zero arguments
+template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *)) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (*f)(Walker);
+  }
+};
+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 ========================== //
+
+// zero arguments
+template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (copy->*f)(Walker);
+  }
+};
+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 ========================== //
+
+// zero arguments
+template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)())
+      (copy->*f)(Walker);
+  }
+};
+template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () const ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)())
+      (copy->*f)(Walker);
+  }
+};
+template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)())
+      (copy->*f)(Walker);
+  }
+};
+template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) () const ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)())
+      (copy->*f)(Walker);
+  }
+};
 // one argument
-template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T), T t )
-{
-  res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)(t);
-  }
-  return result;
-};
-template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T) const, T t )
-{
-  res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)(t);
-  }
-  return result;
-};
-template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T), T t ) const
-{
-  res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)(t);
-  }
-  return result;
-};
-template <typename res, typename typ, typename T> res molecule::SumPerAtom(res (typ::*f)(T) const, T t ) const
-{
-  res result = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    result += (Walker->*f)(t);
-  }
-  return result;
-};
-
-
-// ================== Acting with each Atoms on same molecule ========================== //
+template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T) const, T t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T), T t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T) const, T t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t))
+      (copy->*f)(Walker);
+  }
+};
+// two arguments
+template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u))
+      (copy->*f)(Walker);
+  }
+};
+// three arguments
+template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u,v))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u,v))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u,v))
+      (copy->*f)(Walker);
+  }
+};
+template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if ((Walker->*condition)(t,u,v))
+      (copy->*f)(Walker);
+  }
+};
+
+// ================== Acting on all Atoms ========================== //
 
 // zero arguments
-template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *))
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*f)(Walker);
-  }
-};
-template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *) const)
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*f)(Walker);
-  }
-};
-template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *)) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*f)(Walker);
-  }
-};
-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 ========================== //
-
-// zero arguments
-template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy)
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (copy->*f)(Walker);
-  }
-};
-template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const , molecule *copy)
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (copy->*f)(Walker);
-  }
-};
-template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (copy->*f)(Walker);
-  }
-};
-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 ========================== //
-
-// zero arguments
-template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)())
-      (copy->*f)(Walker);
+template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)()) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)();
+  }
+};
+template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    (Walker->*f)();
   }
 };
 // one argument
-template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t))
-      (copy->*f)(Walker);
-  }
-};
-// two arguments
-template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u))
-      (copy->*f)(Walker);
-  }
-};
-// three arguments
-template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v )
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if ((Walker->*condition)(t,u,v))
-      (copy->*f)(Walker);
-  }
-};
-
-// ================== Acting on all Atoms ========================== //
-
-// zero arguments
-template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)())
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)();
-  }
-};
-template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const)
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)();
-  }
-};
-template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)()) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)();
-  }
-};
-template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)();
-  }
-};
-// one argument
-template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t )
+template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t ) const
 {
   atom *Walker = start;
@@ -398,5 +354,5 @@
   }
 };
-template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t )
+template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t ) const
 {
   atom *Walker = start;
@@ -406,22 +362,6 @@
   }
 };
-template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t ) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t);
-  }
-};
-template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t ) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t);
-  }
-};
 // two argument
-template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*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 ) const
 {
   atom *Walker = start;
@@ -431,5 +371,5 @@
   }
 };
-template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*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 ) const
 {
   atom *Walker = start;
@@ -439,22 +379,6 @@
   }
 };
-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;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u);
-  }
-};
-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;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u);
-  }
-};
 // three argument
-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)
+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;
@@ -464,5 +388,5 @@
   }
 };
-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)
+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;
@@ -472,22 +396,6 @@
   }
 };
-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;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, 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) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v);
-  }
-};
 // four arguments
-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)
+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;
@@ -497,5 +405,5 @@
   }
 };
-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)
+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;
@@ -505,25 +413,9 @@
   }
 };
-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;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v, 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) const
-{
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (Walker->*f)(t, u, v, w);
-  }
-};
 
 // ===================== Accessing arrays indexed by some integer for each atom ======================
 
 // for atom ints
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *) )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *) ) const
 {
   atom *Walker = start;
@@ -534,5 +426,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T value ) const
 {
   atom *Walker = start;
@@ -542,5 +434,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T *value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T *value ) const
 {
   atom *Walker = start;
@@ -551,5 +443,5 @@
 };
 // for element ints
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *) )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *) ) const
 {
   atom *Walker = start;
@@ -560,5 +452,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T value ) const
 {
   atom *Walker = start;
@@ -568,5 +460,5 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T *value )
+template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T *value ) const
 {
   atom *Walker = start;
@@ -577,5 +469,5 @@
 };
 
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(Vector &), Vector atom::*value )
+template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ atom::*value ) const
 {
   atom *Walker = start;
@@ -585,5 +477,13 @@
   }
 };
-template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(Vector &), Vector &vect )
+template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ atom::*value ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    array[(Walker->*index)] = (Walker->*Setor) (Walker->*value);
+  }
+};
+template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ &vect ) const
 {
   atom *Walker = start;
@@ -593,6 +493,13 @@
   }
 };
-
-template <typename T, typename typ, typename typ2> void molecule::SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value )
+template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ &vect ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    array[(Walker->*index)] = (Walker->*Setor) (vect);
+  }
+};
+template <typename T, typename typ, typename typ2> void molecule::SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value ) const
 {
   atom *Walker = start;
@@ -604,5 +511,5 @@
 };
 
-template <typename T, typename typ> void molecule::SetAtomValueToValue ( T value, T typ::*ptr )
+template <typename T, typename typ> void molecule::SetAtomValueToValue ( T value, T typ::*ptr ) const
 {
   atom *Walker = start;
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/vector.cpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -1276,5 +1276,5 @@
  * @param three vectors forming the matrix that defines the shape of the parallelpiped
  */
-bool Vector::IsInParallelepiped(Vector offset, double *parallelepiped)
+bool Vector::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
 {
   Vector a;
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision 831a14fccc3c449e545ce45a4036bbbfeb01b80c)
+++ src/vector.hpp	(revision b453f9633192bf031a859545dcb7c55476ff6a0f)
@@ -77,5 +77,5 @@
   void AskPosition(double *cell_size, bool check);
   bool Output(ofstream *out) const;
-  bool IsInParallelepiped(Vector offset, double *parallelepiped);
+  bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
   void WrapPeriodically(const double *M, const double *Minv);
 };
