Index: src/Analysis/analysis_correlation.cpp
===================================================================
--- src/Analysis/analysis_correlation.cpp	(revision a58c16f47065e9446b7916a407271c6544180022)
+++ src/Analysis/analysis_correlation.cpp	(revision f3b5971ec8cc96538608d12dd2998861f1e187db)
@@ -125,5 +125,5 @@
  * \return range with [min, max]
  */
-range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms)
+range<size_t> getMaximumTrajectoryBounds(const std::vector<const atom *> &atoms)
 {
   // get highest trajectory size
@@ -133,5 +133,5 @@
   size_t max_timesteps = std::numeric_limits<size_t>::min();
   size_t min_timesteps = std::numeric_limits<size_t>::max();
-  BOOST_FOREACH(atom *_atom, atoms) {
+  BOOST_FOREACH(const atom *_atom, atoms) {
     if (_atom->getTrajectorySize() > max_timesteps)
       max_timesteps  = _atom->getTrajectorySize();
@@ -149,17 +149,12 @@
  * \return map with orientation vector for each atomic id given in \a atoms.
  */
-std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<molecule *> &molecules)
+std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<const molecule *> &molecules)
 {
   // get zero orientation for each molecule.
   LOG(0,"STATUS: Calculating dipoles for current time step ...");
   std::map<atomId_t, Vector> ZeroVector;
-  BOOST_FOREACH(molecule *_mol, molecules) {
-    const Vector Dipole =
-        getDipole(
-            const_cast<const molecule *>(_mol)->begin(),
-            const_cast<const molecule *>(_mol)->end());
-    for(molecule::const_iterator iter = const_cast<const molecule *>(_mol)->begin();
-        iter != const_cast<const molecule *>(_mol)->end();
-        ++iter)
+  BOOST_FOREACH(const molecule *_mol, molecules) {
+    const Vector Dipole = getDipole(_mol->begin(),_mol->end());
+    for(molecule::const_iterator iter = _mol->begin(); iter != _mol->end(); ++iter)
       ZeroVector[(*iter)->getId()] = Dipole;
     LOG(2,"INFO: Zero alignment for molecule " << _mol->getId() << " is " << Dipole);
@@ -205,5 +200,6 @@
   World::getInstance().clearMoleculeSelection();
   World::getInstance().selectAllMolecules(MoleculeByFormula(DipoleFormula));
-  std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
+  std::vector<const molecule *> molecules =
+      const_cast<const World &>(World::getInstance()).getSelectedMolecules();
   LOG(1,"INFO: There are " << molecules.size() << " molecules for time step " << timestep << ".");
 
@@ -212,23 +208,20 @@
   size_t i=0;
   size_t Counter_rejections = 0;
-  BOOST_FOREACH(molecule *_mol, molecules) {
-    const Vector Dipole =
-        getDipole(
-            const_cast<const molecule *>(_mol)->begin(),
-            const_cast<const molecule *>(_mol)->end());
+  BOOST_FOREACH(const molecule *_mol, molecules) {
+    const Vector Dipole = getDipole(_mol->begin(),_mol->end());
     LOG(3,"INFO: Dipole vector at time step " << timestep << " for for molecule "
         << _mol->getId() << " is " << Dipole);
     // check that all atoms are valid (zeroVector known)
-    molecule::const_iterator iter = const_cast<const molecule *>(_mol)->begin();
-    for(; iter != const_cast<const molecule *>(_mol)->end(); ++iter) {
+    molecule::const_iterator iter = _mol->begin();
+    for(; iter != _mol->end(); ++iter) {
       if (!ZeroVector.count((*iter)->getId()))
         break;
     }
-    if (iter != const_cast<const molecule *>(_mol)->end()) {
+    if (iter != _mol->end()) {
       ELOG(2, "Skipping molecule " << _mol->getName() << " as not all atoms have a valid zeroVector.");
       ++Counter_rejections;
       continue;
     } else
-      iter = const_cast<const molecule *>(_mol)->begin();
+      iter = _mol->begin();
     std::map<atomId_t, Vector>::const_iterator zeroValue = ZeroVector.find((*iter)->getId()); //due to iter is const
     double angle = 0.;
@@ -274,5 +267,6 @@
  * \return Map of doubles with values the pair of the two atoms.
  */
-DipoleCorrelationMap *DipoleCorrelation(std::vector<molecule *> &molecules)
+DipoleCorrelationMap *DipoleCorrelation(
+    const std::vector<const molecule *> &molecules)
 {
   Info FunctionInfo(__func__);
@@ -286,19 +280,13 @@
   }
 
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin();
+  for (std::vector<const molecule *>::const_iterator MolWalker = molecules.begin();
       MolWalker != molecules.end(); ++MolWalker) {
     LOG(2, "INFO: Current molecule is " << (*MolWalker)->getId() << ".");
-    const Vector Dipole =
-        getDipole(
-            const_cast<const molecule *>(*MolWalker)->begin(),
-            const_cast<const molecule *>(*MolWalker)->end());
-    std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker;
-    for (++MolOtherWalker;
-        MolOtherWalker != molecules.end();
+    const Vector Dipole = getDipole((*MolWalker)->begin(), (*MolWalker)->end());
+    std::vector<const molecule *>::const_iterator MolOtherWalker = MolWalker;
+    for (++MolOtherWalker; MolOtherWalker != molecules.end();
         ++MolOtherWalker) {
       LOG(2, "INFO: Current other molecule is " << (*MolOtherWalker)->getId() << ".");
-      const Vector OtherDipole = getDipole(
-          const_cast<const molecule *>(*MolOtherWalker)->begin(),
-          const_cast<const molecule *>(*MolOtherWalker)->end());
+      const Vector OtherDipole = getDipole((*MolOtherWalker)->begin(), (*MolOtherWalker)->end());
       const double angle = Dipole.Angle(OtherDipole) * (180./M_PI);
       LOG(1, "Angle is " << angle << ".");
@@ -409,5 +397,8 @@
  * \return Map of dobules with values as pairs of atom and the vector
  */
-CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point )
+CorrelationToPointMap *CorrelationToPoint(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Vector *point )
 {
   Info FunctionInfo(__func__);
@@ -421,9 +412,9 @@
   }
 
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+  for (std::vector<const molecule *>::const_iterator MolWalker = molecules.begin();
+      MolWalker != molecules.end();
+      MolWalker++) {
     LOG(2, "Current molecule is " << *MolWalker << ".");
-    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
-        iter != const_cast<const molecule *>(*MolWalker)->end();
-        ++iter) {
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
       LOG(3, "Current atom is " << **iter << ".");
       for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
@@ -453,5 +444,9 @@
  * \return Map of dobules with values as pairs of atom and the vector
  */
-CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] )
+CorrelationToPointMap *PeriodicCorrelationToPoint(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Vector *point,
+    const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -467,11 +462,11 @@
   }
 
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+  for (std::vector<const molecule *>::const_iterator MolWalker = molecules.begin();
+      MolWalker != molecules.end();
+      MolWalker++) {
     RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
     RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     LOG(2, "Current molecule is " << *MolWalker << ".");
-    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
-        iter != const_cast<const molecule *>(*MolWalker)->end();
-        ++iter) {
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
       LOG(3, "Current atom is " << **iter << ".");
       for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
@@ -509,5 +504,9 @@
  * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
  */
-CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC )
+CorrelationToSurfaceMap *CorrelationToSurface(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Tesselation * const Surface,
+    const LinkedCell_deprecated *LC )
 {
   Info FunctionInfo(__func__);
@@ -522,11 +521,11 @@
   }
 
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+  for (std::vector<const molecule *>::const_iterator MolWalker = molecules.begin();
+      MolWalker != molecules.end();
+      MolWalker++) {
     LOG(2, "Current molecule is " << (*MolWalker)->name << ".");
     if ((*MolWalker)->empty())
       LOG(2, "\t is empty.");
-    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
-        iter != const_cast<const molecule *>(*MolWalker)->end();
-        ++iter) {
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
       LOG(3, "\tCurrent atom is " << *(*iter) << ".");
       for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
@@ -563,5 +562,10 @@
  * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
  */
-CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC, const int ranges[NDIM] )
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Tesselation * const Surface,
+    const LinkedCell_deprecated *LC,
+    const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -581,11 +585,11 @@
   double ShortestDistance = 0.;
   BoundaryTriangleSet *ShortestTriangle = NULL;
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+  for (std::vector<const molecule *>::const_iterator MolWalker = molecules.begin();
+      MolWalker != molecules.end();
+      MolWalker++) {
     RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
     RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     LOG(2, "Current molecule is " << *MolWalker << ".");
-    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
-        iter != const_cast<const molecule *>(*MolWalker)->end();
-        ++iter) {
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
       LOG(3, "Current atom is " << **iter << ".");
       for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
Index: src/Analysis/analysis_correlation.hpp
===================================================================
--- src/Analysis/analysis_correlation.hpp	(revision a58c16f47065e9446b7916a407271c6544180022)
+++ src/Analysis/analysis_correlation.hpp	(revision f3b5971ec8cc96538608d12dd2998861f1e187db)
@@ -61,30 +61,73 @@
 /********************************************** declarations *******************************/
 
-range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms);
-std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<molecule *> &molecules);
-
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(const Formula &DipoleFormula, const size_t timestep, const std::map<atomId_t, Vector> &ZeroVector, const enum ResetWorldTime DoTimeReset = DontResetTime);
-DipoleCorrelationMap *DipoleCorrelation(std::vector<molecule *> &molecules);
+range<size_t> getMaximumTrajectoryBounds(
+    const std::vector<const atom *> &atoms);
+std::map<atomId_t, Vector> CalculateZeroAngularDipole(
+    const std::vector<const molecule *> &molecules);
+
+DipoleAngularCorrelationMap *DipoleAngularCorrelation(
+    const Formula &DipoleFormula,
+    const size_t timestep,
+    const std::map<atomId_t, Vector> &ZeroVector,
+    const enum ResetWorldTime DoTimeReset = DontResetTime);
+DipoleCorrelationMap *DipoleCorrelation(
+    const std::vector<const molecule *> &molecules);
 PairCorrelationMap *PairCorrelation(
     const World::ConstAtomComposite &atoms_first,
     const World::ConstAtomComposite &atoms_second,
     const double max_distance);
-CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point );
-CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC );
-CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] );
-CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC, const int ranges[NDIM] );
-int GetBin ( const double value, const double BinWidth, const double BinStart );
-void OutputCorrelation_Header( ofstream * const file );
-void OutputCorrelation_Value( ofstream * const file, BinPairMap::const_iterator &runner );
-void OutputDipoleAngularCorrelation_Header( ofstream * const file );
-void OutputDipoleAngularCorrelation_Value( ofstream * const file, DipoleAngularCorrelationMap::const_iterator &runner );
-void OutputDipoleCorrelation_Header( ofstream * const file );
-void OutputDipoleCorrelation_Value( ofstream * const file, DipoleCorrelationMap::const_iterator &runner );
-void OutputPairCorrelation_Header( ofstream * const file );
-void OutputPairCorrelation_Value( ofstream * const file, PairCorrelationMap::const_iterator &runner );
-void OutputCorrelationToPoint_Header( ofstream * const file );
-void OutputCorrelationToPoint_Value( ofstream * const file, CorrelationToPointMap::const_iterator &runner );
-void OutputCorrelationToSurface_Header( ofstream * const file );
-void OutputCorrelationToSurface_Value( ofstream * const file, CorrelationToSurfaceMap::const_iterator &runner );
+CorrelationToPointMap *CorrelationToPoint(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Vector *point );
+CorrelationToSurfaceMap *CorrelationToSurface(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Tesselation * const Surface,
+    const LinkedCell_deprecated *LC );
+CorrelationToPointMap *PeriodicCorrelationToPoint(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Vector *point, const int ranges[NDIM] );
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface(
+    const std::vector<const molecule *> &molecules,
+    const std::vector<const element *> &elements,
+    const Tesselation * const Surface,
+    const LinkedCell_deprecated *LC,
+    const int ranges[NDIM] );
+int GetBin (
+    const double value,
+    const double BinWidth,
+    const double BinStart );
+void OutputCorrelation_Header(
+    ofstream * const file );
+void OutputCorrelation_Value(
+    ofstream * const file,
+    BinPairMap::const_iterator &runner );
+void OutputDipoleAngularCorrelation_Header(
+    ofstream * const file );
+void OutputDipoleAngularCorrelation_Value(
+    ofstream * const file,
+    DipoleAngularCorrelationMap::const_iterator &runner );
+void OutputDipoleCorrelation_Header(
+    ofstream * const file );
+void OutputDipoleCorrelation_Value(
+    ofstream * const file,
+    DipoleCorrelationMap::const_iterator &runner );
+void OutputPairCorrelation_Header(
+    ofstream * const file );
+void OutputPairCorrelation_Value(
+    ofstream * const file,
+    PairCorrelationMap::const_iterator &runner );
+void OutputCorrelationToPoint_Header(
+    ofstream * const file );
+void OutputCorrelationToPoint_Value(
+    ofstream * const file,
+    CorrelationToPointMap::const_iterator &runner );
+void OutputCorrelationToSurface_Header(
+    ofstream * const file );
+void OutputCorrelationToSurface_Value(
+    ofstream * const file,
+    CorrelationToSurfaceMap::const_iterator &runner );
 
 /** Searches for lowest and highest value in a given map of doubles.
Index: src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision a58c16f47065e9446b7916a407271c6544180022)
+++ src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision f3b5971ec8cc96538608d12dd2998861f1e187db)
@@ -105,5 +105,5 @@
   // init maps
   World::getInstance().selectAllMolecules(AllMolecules());
-  allMolecules = World::getInstance().getSelectedMolecules();
+  allMolecules = const_cast<const World &>(World::getInstance()).getSelectedMolecules();
   CPPUNIT_ASSERT_EQUAL( (size_t) 1, allMolecules.size());
   pointmap = CorrelationToPoint( allMolecules, elements, (const Vector *)point );
Index: src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp	(revision a58c16f47065e9446b7916a407271c6544180022)
+++ src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp	(revision f3b5971ec8cc96538608d12dd2998861f1e187db)
@@ -41,5 +41,5 @@
 private:
 
-      std::vector<molecule *> allMolecules;
+      std::vector<const molecule *> allMolecules;
       molecule *TestMolecule;
       const element *hydrogen;
Index: src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision a58c16f47065e9446b7916a407271c6544180022)
+++ src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision f3b5971ec8cc96538608d12dd2998861f1e187db)
@@ -145,5 +145,5 @@
 
   World::getInstance().selectAllMolecules(AllMolecules());
-  allMolecules = World::getInstance().getSelectedMolecules();
+  allMolecules = const_cast<const World &>(World::getInstance()).getSelectedMolecules();
   CPPUNIT_ASSERT_EQUAL( (size_t) 2, allMolecules.size());
 
Index: src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp	(revision a58c16f47065e9446b7916a407271c6544180022)
+++ src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp	(revision f3b5971ec8cc96538608d12dd2998861f1e187db)
@@ -49,5 +49,5 @@
 private:
 
-      std::vector<molecule *> allMolecules;
+      std::vector<const molecule *> allMolecules;
       molecule *TestSurfaceMolecule;
       const element *hydrogen;
