Index: src/Analysis/analysis_correlation.cpp
===================================================================
--- src/Analysis/analysis_correlation.cpp	(revision 505d059b8f5c10ec0e270a2fbd0cbe8462d545d9)
+++ src/Analysis/analysis_correlation.cpp	(revision e658783525483171aa12ab81bf75d42b28dba211)
@@ -32,4 +32,6 @@
 #include "CodePatterns/Log.hpp"
 #include "CodePatterns/Verbose.hpp"
+#include "Descriptors/AtomOfMoleculeSelectionDescriptor.hpp"
+#include "Descriptors/MoleculeFormulaDescriptor.hpp"
 #include "Descriptors/MoleculeOfAtomSelectionDescriptor.hpp"
 #include "Formula.hpp"
@@ -100,5 +102,5 @@
  * \return range with [min, max]
  */
-range<size_t> getMaximumTrajectoryBounds(std::vector<atom *> &atoms)
+range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms)
 {
   // get highest trajectory size
@@ -121,16 +123,11 @@
 
 /** Calculates the angular dipole zero orientation from current time step.
- * \param atoms vector of atoms to calculate it for
+ * \param molecules vector of molecules to calculate dipoles of
  * \return map with orientation vector for each atomic id given in \a atoms.
  */
-std::map<atomId_t, Vector> CalculateZeroAngularDipole(std::vector<atom *> &atoms)
-{
-  // calculate molecules for this time step
-  std::set<molecule *> molecules;
-  BOOST_FOREACH(atom *_atom, atoms)
-    molecules.insert(_atom->getMolecule());
-
+std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<molecule *> &molecules)
+{
   // get zero orientation for each molecule.
-  LOG(0,"STATUS: Calculating dipoles for first time step ...");
+  LOG(0,"STATUS: Calculating dipoles for current time step ...");
   std::map<atomId_t, Vector> ZeroVector;
   BOOST_FOREACH(molecule *_mol, molecules) {
@@ -154,12 +151,11 @@
  *  \a ZeroVector)
  * \param ZeroVector map with Zero orientation vector for each atom in \a atoms.
- *   Is filled from initial time step if size of map does not match size of \a atoms.
  * \param DontResetTime don't reset time to old value (triggers re-creation of bond system)
  * \return Map of doubles with values the pair of the two atoms.
  */
 DipoleAngularCorrelationMap *DipoleAngularCorrelation(
-    std::vector<atom *> &atoms,
+    const Formula &DipoleFormula,
     const size_t timestep,
-    std::map<atomId_t, Vector> &ZeroVector,
+    const std::map<atomId_t, Vector> &ZeroVector,
     const enum ResetWorldTime DoTimeReset
     )
@@ -167,12 +163,4 @@
   Info FunctionInfo(__func__);
   DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;
-
-  // get zero orientation for each molecule if not given
-  if (ZeroVector.size() != atoms.size()) {
-    LOG(1, "INFO: Mismatch size of zero orientation map ("
-        << ZeroVector.size() << ") and atom vector ("<< atoms.size() << +").");
-    ZeroVector.clear();
-    ZeroVector = CalculateZeroAngularDipole(atoms);
-  }
 
   unsigned int oldtime = 0;
@@ -187,26 +175,37 @@
 
   // get all molecules for this time step
-  LOG(0,"STATUS: Gathering molecules from " << atoms.size() << " atoms for time step " << timestep << " ...");
-  std::set<molecule *> molecules;
-  BOOST_FOREACH(atom *_atom, atoms)
-    molecules.insert(_atom->getMolecule());
+  World::getInstance().clearMoleculeSelection();
+  World::getInstance().selectAllMolecules(MoleculeByFormula(DipoleFormula));
+  World::getInstance().clearAtomSelection();
+  World::getInstance().selectAllAtoms(AtomsByMoleculeSelection());
+  std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
+  LOG(1,"STATUS: Gathering " << molecules.size() << " molecules for time step " << timestep << ".");
 
   // calculate dipoles for each
-  LOG(0,"STATUS: Calculating dipoles for " << molecules.size() << " molecules for time step " << timestep << " ...");
+  LOG(1,"STATUS: Calculating dipoles for time step " << timestep << " ...");
   size_t i=0;
   BOOST_FOREACH(molecule *_mol, molecules) {
     const Vector Dipole = getDipole(_mol->begin(), _mol->end());
-    LOG(2,"INFO: Dipole vector at time step " << timestep << " for molecule "
+    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 = _mol->begin();
-    ASSERT(ZeroVector.count((*iter)->getId()),
-        "DipoleAngularCorrelation() - ZeroVector for atom "+toString(**iter)+" not present.");
+    for(; iter != _mol->end(); ++iter) {
+      if (!ZeroVector.count((*iter)->getId()))
+        break;
+    }
+    if (iter != _mol->end()) {
+      ELOG(2, "Skipping molecule " << _mol->getName() << " as not all atoms have a valid zeroVector.");
+      continue;
+    } else
+      iter = _mol->begin();
+    std::map<atomId_t, Vector>::const_iterator zeroValue = ZeroVector.find((*iter)->getId()); //due to iter is const
     double angle = 0.;
     LOG(2, "INFO: ZeroVector of first atom " << **iter << " is "
-        << ZeroVector[(*iter)->getId()] << ".");
+        << zeroValue->second << ".");
     LOG(4, "INFO: Squared norm of difference vector is "
-        << (ZeroVector[(*iter)->getId()] - Dipole).NormSquared() << ".");
-    if ((ZeroVector[(*iter)->getId()] - Dipole).NormSquared() > MYEPSILON)
-      angle = Dipole.Angle(ZeroVector[(*iter)->getId()]) * (180./M_PI);
+        << (zeroValue->second - Dipole).NormSquared() << ".");
+    if ((zeroValue->second - Dipole).NormSquared() > MYEPSILON)
+      angle = Dipole.Angle(zeroValue->second) * (180./M_PI);
     else
       LOG(2, "INFO: Both vectors (almost) coincide, numerically unstable, angle set to zero.");
@@ -216,5 +215,5 @@
     ++i;
   }
-  LOG(0,"STATUS: Done with calculating dipoles.");
+  LOG(1,"STATUS: Done with calculating dipoles.");
 
   if (DoTimeReset == DoResetTime) {
Index: src/Analysis/analysis_correlation.hpp
===================================================================
--- src/Analysis/analysis_correlation.hpp	(revision 505d059b8f5c10ec0e270a2fbd0cbe8462d545d9)
+++ src/Analysis/analysis_correlation.hpp	(revision e658783525483171aa12ab81bf75d42b28dba211)
@@ -60,8 +60,8 @@
 /********************************************** declarations *******************************/
 
-range<size_t> getMaximumTrajectoryBounds(std::vector<atom *> &atoms);
-std::map<atomId_t, Vector> CalculateZeroAngularDipole(std::vector<atom *> &atoms);
-
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<atom *> &atoms, const size_t timestep, std::map<atomId_t, Vector> &ZeroVector, const enum ResetWorldTime DoTimeReset = DontResetTime);
+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);
 PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements);
