Index: src/Analysis/Makefile.am
===================================================================
--- src/Analysis/Makefile.am	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/Makefile.am	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,52 @@
+# PLEASE adhere to the alphabetical ordering in this Makefile!
+# Also indentation by a single tab
+
+ANALYSISSOURCE = \
+	Analysis/analysis_bonds.cpp \
+	Analysis/analysis_correlation.cpp
+				  
+ANALYSISHEADER = \
+	Analysis/analysis_bonds.hpp \
+	Analysis/analysis_correlation.hpp
+
+
+noinst_LTLIBRARIES += libMolecuilderAnalysis.la
+libMolecuilderAnalysis_la_includedir = $(includedir)/MoleCuilder/Analysis/
+
+nobase_libMolecuilderAnalysis_la_include_HEADERS = ${ANALYSISHEADER}
+
+## Define the source file list for the "libexample-@MOLECUILDER_API_VERSION@.la"
+## target.  Note that @MOLECUILDER_API_VERSION@ is not interpreted by Automake and
+## will therefore be treated as if it were literally part of the target name,
+## and the variable name derived from that.
+## The file extension .cc is recognized by Automake, and makes it produce
+## rules which invoke the C++ compiler to produce a libtool object file (.lo)
+## from each source file.  Note that it is not necessary to list header files
+## which are already listed elsewhere in a _HEADERS variable assignment.
+libMolecuilderAnalysis_la_SOURCES = ${ANALYSISSOURCE}
+
+## Instruct libtool to include ABI version information in the generated shared
+## library file (.so).  The library ABI version is defined in configure.ac, so
+## that all version information is kept in one place.
+libMolecuilderAnalysis_la_LDFLAGS = -version-info $(MOLECUILDER_SO_VERSION)
+
+## The generated configuration header is installed in its own subdirectory of
+## $(libdir).  The reason for this is that the configuration information put
+## into this header file describes the target platform the installed library
+## has been built for.  Thus the file must not be installed into a location
+## intended for architecture-independent files, as defined by the Filesystem
+## Hierarchy Standard (FHS).
+## The nodist_ prefix instructs Automake to not generate rules for including
+## the listed files in the distribution on 'make dist'.  Files that are listed
+## in _HEADERS variables are normally included in the distribution, but the
+## configuration header file is generated at configure time and should not be
+## shipped with the source tarball.
+#libMolecuilderAnalysis_libincludedir = $(libdir)/MoleCuilder/include
+#nodist_libMolecuilderAnalysis_libinclude_HEADERS = $(top_builddir)/libmolecuilder_config.h
+
+## Install the generated pkg-config file (.pc) into the expected location for
+## architecture-dependent package configuration information.  Occasionally,
+## pkg-config files are also used for architecture-independent data packages,
+## in which case the correct install location would be $(datadir)/pkgconfig.
+#pkgconfigdir = $(libdir)/pkgconfig
+#pkgconfig_DATA = $(top_builddir)/MoleCuilder.pc
Index: src/Analysis/analysis_bonds.cpp
===================================================================
--- src/Analysis/analysis_bonds.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/analysis_bonds.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,299 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * analysis_bonds.cpp
+ *
+ *  Created on: Nov 7, 2009
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "analysis_bonds.hpp"
+#include "atom.hpp"
+#include "Bond/bond.hpp"
+#include "element.hpp"
+#include "CodePatterns/Info.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "CodePatterns/Log.hpp"
+#include "molecule.hpp"
+
+/** Calculates the min, mean and maximum bond counts for the given molecule.
+ * \param *mol molecule with atoms and atom::ListOfBonds
+ * \param &Min minimum count on return
+ * \param &Mean mean count on return
+ * \param &Max maximum count on return
+ */
+void GetMaxMinMeanBondCount(const molecule * const mol, double &Min, double &Mean, double &Max)
+{
+  Min = 2e+6;
+  Max = -2e+5;
+  Mean = 0.;
+
+  int AtomCount = 0;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    const BondList& ListOfBonds = (*iter)->getListOfBonds();
+    const int count = ListOfBonds.size();
+    if (Max < count)
+      Max = count;
+    if (Min > count)
+      Min = count;
+    Mean += count;
+    AtomCount++;
+  }
+  if (((int)Mean % 2) != 0)
+    DoeLog(1) && (eLog()<< Verbose(1) << "Something is wrong with the bond structure, the number of bonds is not even!" << endl);
+  Mean /= (double)AtomCount;
+};
+
+/** Calculates the min and max bond distance of all atoms of two given elements.
+ * \param *mol molecule with atoms
+ * \param *type1 one element
+ * \param *type2 other element
+ * \param &Min minimum distance on return, 0 if no bond between the two elements
+ * \param &Mean mean distance (i.e. sum of distance for matching element pairs, divided by number) on return, 0 if no bond between the two elements
+ * \param &Max maximum distance on return, 0 if no bond between the two elements
+ */
+void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, const element *type1, const element *type2, double &Min, double &Mean, double &Max)
+{
+  Min = 2e+6;
+  Mean = 0.;
+  Max = -2e+6;
+
+  int AtomNo = 0;
+  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+    if ((*iter)->getType() == type1) {
+      const BondList& ListOfBonds = (*iter)->getListOfBonds();
+      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+          BondRunner != ListOfBonds.end();
+          BondRunner++)
+        if ((*BondRunner)->GetOtherAtom((*iter))->getType() == type2) {
+          const double distance = (*BondRunner)->GetDistanceSquared();
+          if (Min > distance)
+            Min = distance;
+          if (Max < distance)
+            Max = distance;
+          Mean += sqrt(distance);
+          AtomNo++;
+        }
+    }
+  }
+  if (Max < 0) {
+    Max = Min = 0.;
+  } else {
+    Max = sqrt(Max);
+    Min = sqrt(Min);
+    Mean = Mean/(double)AtomNo;
+  }
+};
+
+/** Calculate the angle between \a *first and \a *origin and \a *second and \a *origin.
+ * \param *first first Vector
+ * \param *origin origin of angle taking
+ * \param *second second Vector
+ * \return angle between \a *first and \a *second, both relative to origin at \a *origin.
+ */
+double CalculateAngle(const Vector &first, const Vector &central, const Vector &second)
+{
+  Vector OHBond;
+  Vector OOBond;
+
+  OHBond = first - central;
+  OOBond = second - central;
+  const double angle = OHBond.Angle(OOBond);
+  return angle;
+};
+
+/** Checks whether the angle between \a *Oxygen and \a *Hydrogen and \a *Oxygen and \a *OtherOxygen is less than 30 degrees.
+ * Note that distance criterion is not checked.
+ * \param *Oxygen first oxygen atom, bonded to \a *Hydrogen
+ * \param *Hydrogen hydrogen bonded to \a *Oxygen
+ * \param *OtherOxygen other oxygen atom
+ * \return true - angle criteria fulfilled, false - criteria not fulfilled, angle greater than 30 degrees.
+ */
+bool CheckHydrogenBridgeBondAngle(atom *Oxygen, atom *Hydrogen, atom *OtherOxygen)
+{
+  Info FunctionInfo(__func__);
+
+  // check angle
+  if (CalculateAngle(Hydrogen->getPosition(), Oxygen->getPosition(), OtherOxygen->getPosition()) < M_PI*(30./180.)) {
+    return true;
+  } else {
+    return false;
+  }
+};
+
+/** Counts the number of hydrogen bridge bonds.
+ * With \a *InterfaceElement an extra element can be specified that identifies some boundary.
+ * Then, counting is for the h-bridges that connect to interface only.
+ * \param *molecules molecules to count bonds
+ * \param *InterfaceElement or NULL
+ * \param *Interface2Element or NULL
+ */
+int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL, const element * Interface2Element = NULL)
+{
+  int count = 0;
+  int OtherHydrogens = 0;
+  double Otherangle = 0.;
+  bool InterfaceFlag = false;
+  bool Interface2Flag = false;
+  bool OtherHydrogenFlag = true;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); ++MolWalker) {
+    molecule::iterator Walker = (*MolWalker)->begin();
+    for(;Walker!=(*MolWalker)->end();++Walker){
+      for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); ++MolRunner) {
+        molecule::iterator Runner = (*MolRunner)->begin();
+        for(;Runner!=(*MolRunner)->end();++Runner){
+          if (((*Walker)->getType()->getAtomicNumber()  == 8) && ((*Runner)->getType()->getAtomicNumber()  == 8)) {
+            // check distance
+            const double distance = (*Runner)->DistanceSquared(*(*Walker));
+            if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means  different atoms
+              // on other atom(Runner) we check for bond to interface element and
+              // check that O-O line is not in between the shanks of the two connected hydrogens (Otherangle > 104.5)
+              OtherHydrogenFlag = true;
+              Otherangle = 0.;
+              OtherHydrogens = 0;
+              InterfaceFlag = (InterfaceElement == NULL);
+              Interface2Flag = (Interface2Element == NULL);
+              const BondList& ListOfBonds = (*Runner)->getListOfBonds();
+              for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+                  BondRunner != ListOfBonds.end();
+                  BondRunner++) {
+                atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
+                // if hydrogen, check angle to be greater(!) than 30 degrees
+                if (OtherAtom->getType()->getAtomicNumber() == 1) {
+                  const double angle = CalculateAngle(OtherAtom->getPosition(), (*Runner)->getPosition(), (*Walker)->getPosition());
+                  OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
+                  Otherangle += angle;
+                  OtherHydrogens++;
+                }
+                InterfaceFlag = InterfaceFlag || (OtherAtom->getType() == InterfaceElement);
+                Interface2Flag = Interface2Flag || (OtherAtom->getType() == Interface2Element);
+              }
+              DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl);
+              switch (OtherHydrogens) {
+                case 0:
+                case 1:
+                  break;
+                case 2:
+                  OtherHydrogenFlag = OtherHydrogenFlag && (Otherangle > M_PI*(104.5/180.) + MYEPSILON);
+                  break;
+                default: // 3 or more hydrogens ...
+                  OtherHydrogenFlag = false;
+                  break;
+              }
+              if (InterfaceFlag && Interface2Flag && OtherHydrogenFlag) {
+                // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
+                const BondList& ListOfBonds = (*Walker)->getListOfBonds();
+                for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+                    BondRunner != ListOfBonds.end();
+                    BondRunner++) {
+                  atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
+                  if (OtherAtom->getType()->getAtomicNumber() == 1) {
+                    // check angle
+                    if (CheckHydrogenBridgeBondAngle(*Walker, OtherAtom, *Runner)) {
+                      DoLog(1) && (Log() << Verbose(1) << (*Walker)->getName() << ", " << OtherAtom->getName() << " and " << (*Runner)->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(OtherAtom->getPosition(), (*Walker)->getPosition(), (*Runner)->getPosition())*(180./M_PI) << "." << endl);
+                      count++;
+                      break;
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  return count;
+}
+
+/** Counts the number of bonds between two given elements.
+ * \param *molecules list of molecules with all atoms
+ * \param *first pointer to first element
+ * \param *second pointer to second element
+ * \return number of found bonds (\a *first-\a *second)
+ */
+int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second)
+{
+  int count = 0;
+
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
+    molecule::iterator Walker = (*MolWalker)->begin();
+    for(;Walker!=(*MolWalker)->end();++Walker){
+      atom * theAtom = *Walker;
+      if ((theAtom->getType() == first) || (theAtom->getType() == second)) {  // first element matches
+        const BondList& ListOfBonds = theAtom->getListOfBonds();
+        for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+            BondRunner != ListOfBonds.end();
+            BondRunner++) {
+          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
+          if (((OtherAtom->getType() == first) || (OtherAtom->getType() == second)) && (theAtom->getNr() < OtherAtom->getNr())) {
+            count++;
+            DoLog(1) && (Log() << Verbose(1) << *first << "-" << *second << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl);
+          }
+        }
+      }
+    }
+  }
+  return count;
+};
+
+/** Counts the number of bonds between three given elements.
+ * Note that we do not look for arbitrary sequence of given bonds, but \a *second will be the central atom and we check
+ * whether it has bonds to both \a *first and \a *third.
+ * \param *molecules list of molecules with all atoms
+ * \param *first pointer to first element
+ * \param *second pointer to second element
+ * \param *third pointer to third element
+ * \return number of found bonds (\a *first-\a *second-\a *third, \a *third-\a *second-\a *first, respectively)
+ */
+int CountBondsOfThree(MoleculeListClass * const molecules, const element * const first, const element * const second, const element * const third)
+{
+  int count = 0;
+  bool MatchFlag[2];
+  bool result = false;
+  const element * ElementArray[2];
+  ElementArray[0] = first;
+  ElementArray[1] = third;
+
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
+    molecule::iterator Walker = (*MolWalker)->begin();
+    for(;Walker!=(*MolWalker)->end();++Walker){
+      atom *theAtom = *Walker;
+      if (theAtom->getType() == second) {  // first element matches
+        for (int i=0;i<2;i++)
+          MatchFlag[i] = false;
+        const BondList& ListOfBonds = theAtom->getListOfBonds();
+        for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+            BondRunner != ListOfBonds.end();
+            BondRunner++) {
+          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
+          for (int i=0;i<2;i++)
+            if ((!MatchFlag[i]) && (OtherAtom->getType() == ElementArray[i])) {
+              MatchFlag[i] = true;
+              break;  // each bonding atom can match at most one element we are looking for
+            }
+        }
+        result = true;
+        for (int i=0;i<2;i++) // gather results
+          result = result && MatchFlag[i];
+        if (result) { // check results
+          count++;
+          DoLog(1) && (Log() << Verbose(1) << *first << "-" << *second << "-" << *third << " bond found at " << *Walker << "." << endl);
+        }
+      }
+    }
+  }
+  return count;
+};
Index: src/Analysis/analysis_bonds.hpp
===================================================================
--- src/Analysis/analysis_bonds.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/analysis_bonds.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,39 @@
+/*
+ * analysis_bonds.hpp
+ *
+ *  Created on: Nov 7, 2009
+ *      Author: heber
+ */
+
+#ifndef ANALYSIS_BONDS_HPP_
+#define ANALYSIS_BONDS_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/*********************************************** defines ***********************************/
+
+const double HBRIDGEDISTANCE=3.5;   //!< HBridge distance from PCCP Vol 10. 4802-4813
+
+/****************************************** forward declarations *****************************/
+
+class element;
+class MoleculeListClass;
+class molecule;
+
+/********************************************** declarations *******************************/
+
+void GetMaxMinMeanBondCount(const molecule * const mol, double &Min, double &Mean, double &Max);
+void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, const element *type1, const element *type2, double &Min, double &Mean, double &Max);
+
+int CountHydrogenBridgeBonds(MoleculeListClass * const molecules, const element * InterfaceElement, const element * Interface2Element);
+int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second);
+int CountBondsOfThree(MoleculeListClass * const molecules, const element * const first, const element * const second, const element * const third);
+
+#endif /* ANALYSIS_BONDS_HPP_ */
Index: src/Analysis/analysis_correlation.cpp
===================================================================
--- src/Analysis/analysis_correlation.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/analysis_correlation.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,577 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * analysis.cpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include <iostream>
+#include <iomanip>
+
+#include "atom.hpp"
+#include "Bond/bond.hpp"
+#include "Tesselation/BoundaryTriangleSet.hpp"
+#include "Box.hpp"
+#include "element.hpp"
+#include "CodePatterns/Info.hpp"
+#include "CodePatterns/Log.hpp"
+#include "Formula.hpp"
+#include "molecule.hpp"
+#include "Tesselation/tesselation.hpp"
+#include "Tesselation/tesselationhelpers.hpp"
+#include "Tesselation/triangleintersectionlist.hpp"
+#include "World.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "World.hpp"
+#include "Box.hpp"
+
+#include "analysis_correlation.hpp"
+
+/** Calculates the dipole vector of a given atomSet.
+ *
+ *  Note that we use the following procedure as rule of thumb:
+ *   -# go through every bond of the atom
+ *   -# calculate the difference of electronegativities \f$\Delta\mathrm{EN}\f$
+ *   -# if \f$\Delta\mathrm{EN} > 0.5\f$, we align the bond vector in direction of the more negative element
+ *   -# sum up all vectors
+ *   -# finally, divide by the number of summed vectors
+ *
+ * @param atomsbegin begin iterator of atomSet
+ * @param atomsend end iterator of atomset
+ * @return dipole vector
+ */
+Vector getDipole(molecule::const_iterator atomsbegin, molecule::const_iterator atomsend)
+{
+  Vector DipoleVector;
+  size_t SumOfVectors = 0;
+  // go through all atoms
+  for (molecule::const_iterator atomiter = atomsbegin;
+      atomiter != atomsend;
+      ++atomiter) {
+    // go through all bonds
+    const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
+    for (BondList::const_iterator bonditer = ListOfBonds.begin();
+        bonditer != ListOfBonds.end();
+        ++bonditer) {
+      const atom * Otheratom = (*bonditer)->GetOtherAtom(*atomiter);
+      if (Otheratom->getId() > (*atomiter)->getId()) {
+        const double DeltaEN = (*atomiter)->getType()->getElectronegativity()
+            -Otheratom->getType()->getElectronegativity();
+        Vector BondDipoleVector = (*atomiter)->getPosition() - Otheratom->getPosition();
+        // DeltaEN is always positive, gives correct orientation of vector
+        BondDipoleVector.Normalize();
+        BondDipoleVector *= DeltaEN;
+        DipoleVector += BondDipoleVector;
+        SumOfVectors++;
+      }
+    }
+  }
+  DipoleVector *= 1./(double)SumOfVectors;
+  DoLog(1) && (Log() << Verbose(1) << "Resulting dipole vector is " << DipoleVector << std::endl);
+
+  return DipoleVector;
+};
+
+/** Calculates the dipole angular correlation for given molecule type.
+ * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
+ * Angles are given in degrees.
+ * \param *molecules vector of molecules
+ * \return Map of doubles with values the pair of the two atoms.
+ */
+DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules)
+{
+  Info FunctionInfo(__func__);
+  DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;
+//  double distance = 0.;
+//  Box &domain = World::getInstance().getDomain();
+//
+  if (molecules.empty()) {
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
+    return outmap;
+  }
+
+  outmap = new DipoleAngularCorrelationMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin();
+      MolWalker != molecules.end(); ++MolWalker) {
+    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is "
+        << (*MolWalker)->getId() << "." << endl);
+    const Vector Dipole = getDipole((*MolWalker)->begin(), (*MolWalker)->end());
+    std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker;
+    for (++MolOtherWalker;
+        MolOtherWalker != molecules.end();
+        ++MolOtherWalker) {
+      DoLog(2) && (Log() << Verbose(2) << "Current other molecule is "
+          << (*MolOtherWalker)->getId() << "." << endl);
+      const Vector OtherDipole = getDipole((*MolOtherWalker)->begin(), (*MolOtherWalker)->end());
+      const double angle = Dipole.Angle(OtherDipole) * (180./M_PI);
+      DoLog(1) && (Log() << Verbose(1) << "Angle is " << angle << "." << endl);
+      outmap->insert ( make_pair (angle, make_pair ((*MolWalker), (*MolOtherWalker)) ) );
+    }
+  }
+  return outmap;
+};
+
+
+/** Calculates the pair correlation between given elements.
+ * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
+ * \param *molecules vector of molecules
+ * \param &elements vector of elements to correlate
+ * \return Map of doubles with values the pair of the two atoms.
+ */
+PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements)
+{
+  Info FunctionInfo(__func__);
+  PairCorrelationMap *outmap = new PairCorrelationMap;
+  double distance = 0.;
+  Box &domain = World::getInstance().getDomain();
+
+  if (molecules.empty()) {
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
+    return outmap;
+  }
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
+    (*MolWalker)->doCountAtoms();
+
+  // create all possible pairs of elements
+  set <pair<const element *,const element *> > PairsOfElements;
+  if (elements.size() >= 2) {
+    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
+      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
+        if (type1 != type2) {
+          PairsOfElements.insert( make_pair(*type1,*type2) );
+          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
+        }
+  } else if (elements.size() == 1) { // one to all are valid
+    const element *elemental = *elements.begin();
+    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
+    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
+  } else { // all elements valid
+    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
+  }
+
+  outmap = new PairCorrelationMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
+    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
+      for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
+        DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+        for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
+          DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
+          if ((*iter)->getId() < (*runner)->getId()){
+            for (set <pair<const element *, const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
+              if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
+                distance = domain.periodicDistance((*iter)->getPosition(),(*runner)->getPosition());
+                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
+                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
+              }
+          }
+        }
+      }
+    }
+  }
+  return outmap;
+};
+
+/** Calculates the pair correlation between given elements.
+ * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate
+ * \param ranges[NDIM] interval boundaries for the periodic images to scan also
+ * \return Map of doubles with values the pair of the two atoms.
+ */
+PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] )
+{
+  Info FunctionInfo(__func__);
+  PairCorrelationMap *outmap = new PairCorrelationMap;
+  double distance = 0.;
+  int n[NDIM];
+  Vector checkX;
+  Vector periodicX;
+  int Othern[NDIM];
+  Vector checkOtherX;
+  Vector periodicOtherX;
+
+  if (molecules.empty()) {
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
+    return outmap;
+  }
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
+    (*MolWalker)->doCountAtoms();
+
+  // create all possible pairs of elements
+  set <pair<const element *,const element *> > PairsOfElements;
+  if (elements.size() >= 2) {
+    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
+      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
+        if (type1 != type2) {
+          PairsOfElements.insert( make_pair(*type1,*type2) );
+          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
+        }
+  } else if (elements.size() == 1) { // one to all are valid
+    const element *elemental = *elements.begin();
+    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
+    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
+  } else { // all elements valid
+    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
+  }
+
+  outmap = new PairCorrelationMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
+    RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
+    RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
+    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
+      periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
+      // go through every range in xyz and get distance
+      for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
+        for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
+          for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+            checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
+            for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
+                DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+                for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
+                  DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
+                  if ((*iter)->getId() < (*runner)->getId()){
+                    for (set <pair<const element *,const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
+                      if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
+                        periodicOtherX = FullInverseMatrix * ((*runner)->getPosition()); // x now in [0,1)^3
+                        // go through every range in xyz and get distance
+                        for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
+                          for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
+                            for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
+                              checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
+                              distance = checkX.distance(checkOtherX);
+                              //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
+                              outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
+                            }
+                      }
+                    }
+                  }
+                }
+      }
+    }
+  }
+
+  return outmap;
+};
+
+/** Calculates the distance (pair) correlation between a given element and a point.
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate with point
+ * \param *point vector to the correlation point
+ * \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 )
+{
+  Info FunctionInfo(__func__);
+  CorrelationToPointMap *outmap = new CorrelationToPointMap;
+  double distance = 0.;
+  Box &domain = World::getInstance().getDomain();
+
+  if (molecules.empty()) {
+    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
+    return outmap;
+  }
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
+    (*MolWalker)->doCountAtoms();
+  outmap = new CorrelationToPointMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
+      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          distance = domain.periodicDistance((*iter)->getPosition(),*point);
+          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
+          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
+        }
+    }
+  }
+
+  return outmap;
+};
+
+/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate to point
+ * \param *point vector to the correlation point
+ * \param ranges[NDIM] interval boundaries for the periodic images to scan also
+ * \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] )
+{
+  Info FunctionInfo(__func__);
+  CorrelationToPointMap *outmap = new CorrelationToPointMap;
+  double distance = 0.;
+  int n[NDIM];
+  Vector periodicX;
+  Vector checkX;
+
+  if (molecules.empty()) {
+    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
+    return outmap;
+  }
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
+    (*MolWalker)->doCountAtoms();
+  outmap = new CorrelationToPointMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+    RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
+    RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
+    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
+      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
+          // go through every range in xyz and get distance
+          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
+            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
+              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+                checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
+                distance = checkX.distance(*point);
+                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
+                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
+              }
+        }
+    }
+  }
+
+  return outmap;
+};
+
+/** Calculates the distance (pair) correlation between a given element and a surface.
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate to surface
+ * \param *Surface pointer to Tesselation class surface
+ * \param *LC LinkedCell structure to quickly find neighbouring atoms
+ * \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 *LC )
+{
+  Info FunctionInfo(__func__);
+  CorrelationToSurfaceMap *outmap = new CorrelationToSurfaceMap;
+  double distance = 0;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector centroid;
+
+  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
+    return outmap;
+  }
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
+    (*MolWalker)->doCountAtoms();
+  outmap = new CorrelationToSurfaceMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << (*MolWalker)->name << "." << endl);
+    if ((*MolWalker)->empty())
+      DoLog(2) && (2) && (Log() << Verbose(2) << "\t is empty." << endl);
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      DoLog(3) && (Log() << Verbose(3) << "\tCurrent atom is " << *(*iter) << "." << endl);
+      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          TriangleIntersectionList Intersections((*iter)->getPosition(),Surface,LC);
+          distance = Intersections.GetSmallestDistance();
+          triangle = Intersections.GetClosestTriangle();
+          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
+        }
+    }
+  }
+
+  return outmap;
+};
+
+/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
+ * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
+ * I.e. We multiply the atom::node with the inverse of the domain matrix, i.e. transform it to \f$[0,0^3\f$, then add per
+ * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
+ * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate to surface
+ * \param *Surface pointer to Tesselation class surface
+ * \param *LC LinkedCell structure to quickly find neighbouring atoms
+ * \param ranges[NDIM] interval boundaries for the periodic images to scan also
+ * \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 *LC, const int ranges[NDIM] )
+{
+  Info FunctionInfo(__func__);
+  CorrelationToSurfaceMap *outmap = new CorrelationToSurfaceMap;
+  double distance = 0;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector centroid;
+  int n[NDIM];
+  Vector periodicX;
+  Vector checkX;
+
+  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
+    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
+    return outmap;
+  }
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
+    (*MolWalker)->doCountAtoms();
+  outmap = new CorrelationToSurfaceMap;
+  double ShortestDistance = 0.;
+  BoundaryTriangleSet *ShortestTriangle = NULL;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
+    RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
+    RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
+    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
+      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
+      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+        if ((*type == NULL) || ((*iter)->getType() == *type)) {
+          periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
+          // go through every range in xyz and get distance
+          ShortestDistance = -1.;
+          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
+            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
+              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+                checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
+                TriangleIntersectionList Intersections(checkX,Surface,LC);
+                distance = Intersections.GetSmallestDistance();
+                triangle = Intersections.GetClosestTriangle();
+                if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
+                  ShortestDistance = distance;
+                  ShortestTriangle = triangle;
+                }
+              }
+          // insert
+          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
+          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
+        }
+    }
+  }
+
+  return outmap;
+};
+
+/** Returns the index of the bin for a given value.
+ * \param value value whose bin to look for
+ * \param BinWidth width of bin
+ * \param BinStart first bin
+ */
+int GetBin ( const double value, const double BinWidth, const double BinStart )
+{
+  //Info FunctionInfo(__func__);
+  int bin =(int) (floor((value - BinStart)/BinWidth));
+  return (bin);
+};
+
+
+/** Adds header part that is unique to BinPairMap.
+ *
+ * @param file stream to print to
+ */
+void OutputCorrelation_Header( ofstream * const file )
+{
+  *file << "\tCount";
+};
+
+/** Prints values stored in BinPairMap iterator.
+ *
+ * @param file stream to print to
+ * @param runner iterator pointing at values to print
+ */
+void OutputCorrelation_Value( ofstream * const file, BinPairMap::const_iterator &runner )
+{
+  *file << runner->second;
+};
+
+
+/** Adds header part that is unique to DipoleAngularCorrelationMap.
+ *
+ * @param file stream to print to
+ */
+void OutputDipoleAngularCorrelation_Header( ofstream * const file )
+{
+  *file << "\tAtom1\tAtom2";
+};
+
+/** Prints values stored in DipoleAngularCorrelationMap iterator.
+ *
+ * @param file stream to print to
+ * @param runner iterator pointing at values to print
+ */
+void OutputDipoleAngularCorrelation_Value( ofstream * const file, DipoleAngularCorrelationMap::const_iterator &runner )
+{
+  *file << runner->second.first->getId() << "\t" << runner->second.second->getId();
+};
+
+
+/** Adds header part that is unique to PairCorrelationMap.
+ *
+ * @param file stream to print to
+ */
+void OutputPairCorrelation_Header( ofstream * const file )
+{
+  *file << "\tAtom1\tAtom2";
+};
+
+/** Prints values stored in PairCorrelationMap iterator.
+ *
+ * @param file stream to print to
+ * @param runner iterator pointing at values to print
+ */
+void OutputPairCorrelation_Value( ofstream * const file, PairCorrelationMap::const_iterator &runner )
+{
+  *file << *(runner->second.first) << "\t" << *(runner->second.second);
+};
+
+
+/** Adds header part that is unique to CorrelationToPointMap.
+ *
+ * @param file stream to print to
+ */
+void OutputCorrelationToPoint_Header( ofstream * const file )
+{
+  *file << "\tAtom::x[i]-point.x[i]";
+};
+
+/** Prints values stored in CorrelationToPointMap iterator.
+ *
+ * @param file stream to print to
+ * @param runner iterator pointing at values to print
+ */
+void OutputCorrelationToPoint_Value( ofstream * const file, CorrelationToPointMap::const_iterator &runner )
+{
+  for (int i=0;i<NDIM;i++)
+    *file << "\t" << setprecision(8) << (runner->second.first->at(i) - runner->second.second->at(i));
+};
+
+
+/** Adds header part that is unique to CorrelationToSurfaceMap.
+ *
+ * @param file stream to print to
+ */
+void OutputCorrelationToSurface_Header( ofstream * const file )
+{
+  *file << "\tTriangle";
+};
+
+/** Prints values stored in CorrelationToSurfaceMap iterator.
+ *
+ * @param file stream to print to
+ * @param runner iterator pointing at values to print
+ */
+void OutputCorrelationToSurface_Value( ofstream * const file, CorrelationToSurfaceMap::const_iterator &runner )
+{
+  *file << *(runner->second.first) << "\t" << *(runner->second.second);
+};
Index: src/Analysis/analysis_correlation.hpp
===================================================================
--- src/Analysis/analysis_correlation.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/analysis_correlation.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,191 @@
+/*
+ * analysis.hpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+#ifndef ANALYSIS_HPP_
+#define ANALYSIS_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cmath>
+#include <iomanip>
+#include <fstream>
+
+// STL headers
+#include <map>
+#include <vector>
+
+#include "Helpers/defs.hpp"
+
+#include "atom.hpp"
+#include "CodePatterns/Info.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Helpers/helpers.hpp"
+
+/****************************************** forward declarations *****************************/
+
+class BoundaryTriangleSet;
+class Formula;
+class element;
+class LinkedCell;
+class Tesselation;
+class Vector;
+
+/********************************************** definitions *********************************/
+
+typedef multimap<double, pair<atom *, atom *> > PairCorrelationMap;
+typedef multimap<double, pair<molecule *, molecule *> > DipoleAngularCorrelationMap;
+typedef multimap<double, pair<atom *, const Vector *> > CorrelationToPointMap;
+typedef multimap<double, pair<atom *, BoundaryTriangleSet *> > CorrelationToSurfaceMap;
+typedef map<double, int> BinPairMap;
+
+/********************************************** declarations *******************************/
+
+DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules);
+PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements);
+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 *LC );
+PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] );
+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 *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 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.
+ * \param *map map of doubles to scan
+ * \param &min minimum on return
+ * \param &max maximum on return
+ */
+template <typename T> void GetMinMax ( T *map, double &min, double &max)
+{
+  min = 0.;
+  max = 0.;
+  bool FirstMinFound = false;
+  bool FirstMaxFound = false;
+
+  if (map == NULL) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to min/max, map is NULL!" << endl);
+    performCriticalExit();
+    return;
+  }
+
+  for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
+    if ((min > runner->first) || (!FirstMinFound)) {
+      min = runner->first;
+      FirstMinFound = true;
+    }
+    if ((max < runner->first) || (!FirstMaxFound)) {
+      max = runner->first;
+      FirstMaxFound = true;
+    }
+  }
+};
+
+/** Puts given correlation data into bins of given size (histogramming).
+ * Note that BinStart and BinEnd are desired to fill the complete range, even where Bins are zero. If this is
+ * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater. If a
+ * certain start value is desired, give BinStart and a BinEnd that is smaller than BinStart, then only BinEnd will be
+ * calculated automatically, i.e. BinStart = 0. and BinEnd = -1. .
+ * Also note that the range is given inclusive, i.e. [ BinStart, BinEnd ].
+ * \param *map map of doubles to count
+ * \param BinWidth desired width of the binds
+ * \param BinStart first bin
+ * \param BinEnd last bin
+ * \return Map of doubles (the bins) with counts as values
+ */
+template <typename T> BinPairMap *BinData ( T *map, const double BinWidth, const double BinStart, const double BinEnd )
+{
+  BinPairMap *outmap = new BinPairMap;
+  int bin = 0;
+  double start = 0.;
+  double end = 0.;
+  pair <BinPairMap::iterator, bool > BinPairMapInserter;
+
+  if (map == NULL) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to bin, is NULL!" << endl);
+    performCriticalExit();
+    return outmap;
+  }
+
+  if (BinStart == BinEnd) { // if same, find range ourselves
+    GetMinMax( map, start, end);
+  } else if (BinEnd < BinStart) { // if BinEnd smaller, just look for new max
+    GetMinMax( map, start, end);
+    start = BinStart;
+  } else { // else take both values
+    start = BinStart;
+    end = BinEnd;
+  }
+  for (int runner = 0; runner <= ceil((end-start)/BinWidth); runner++)
+    outmap->insert( pair<double, int> ((double)runner*BinWidth+start, 0) );
+
+  for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
+    bin = GetBin (runner->first, BinWidth, start);
+    BinPairMapInserter = outmap->insert ( pair<double, int> ((double)bin*BinWidth+start, 1) );
+    if (!BinPairMapInserter.second) {  // bin already present, increase
+      BinPairMapInserter.first->second += 1;
+    }
+  }
+
+  return outmap;
+};
+
+/** Prints correlation map of template type to file.
+ * \param *file file to write to
+ * \param *map map to write
+ * \param (*header) pointer to function that adds specific part to header
+ * \param (*value) pointer to function that prints value from given iterator
+ */
+template <typename T>
+void OutputCorrelationMap(
+    ofstream * const file,
+    const T * const map,
+    void (*header)( ofstream * const file),
+    void (*value)( ofstream * const file, typename T::const_iterator &runner)
+    )
+{
+  Info FunctionInfo(__func__);
+  *file << "BinStart\tBinCenter\tBinEnd";
+  header(file);
+  *file << endl;
+  typename T::const_iterator advancer = map->begin();
+  typename T::const_iterator runner = map->begin();
+  if (advancer != map->end())
+    advancer++;
+  for ( ;(advancer != map->end()) && (runner != map->end()); ++advancer, ++runner) {
+    *file << setprecision(8)
+        << runner->first << "\t"
+        << (advancer->first + runner->first)/2. << "\t"
+        << advancer->first << "\t";
+    value(file, runner);
+    *file << endl;
+  }
+  if(runner != map->end()) {
+    *file << setprecision(8) << runner->first << "\t0\t0\t";
+    value(file, runner);
+    *file << endl;
+  }
+};
+
+
+#endif /* ANALYSIS_HPP_ */
Index: src/Analysis/unittests/AnalysisBondsUnitTest.cpp
===================================================================
--- src/Analysis/unittests/AnalysisBondsUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisBondsUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,171 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * AnalysisBondsUnitTest.cpp
+ *
+ *  Created on: Nov 7, 2009
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <iostream>
+#include <stdio.h>
+#include <cstring>
+
+#include "Analysis/analysis_bonds.hpp"
+#include "atom.hpp"
+#include "Bond/bond.hpp"
+#include "element.hpp"
+#include "Graph/BondGraph.hpp"
+#include "molecule.hpp"
+#include "periodentafel.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "World.hpp"
+
+#include "AnalysisBondsUnitTest.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( AnalysisBondsTest );
+
+
+void AnalysisBondsTest::setUp()
+{
+  atom *Walker = NULL;
+
+  // get elements
+  hydrogen = World::getInstance().getPeriode()->FindElement(1);
+  carbon = World::getInstance().getPeriode()->FindElement(6);
+  CPPUNIT_ASSERT(hydrogen != NULL && "could not find element hydrogen");
+  CPPUNIT_ASSERT(carbon != NULL && "could not find element carbon");
+
+  // construct molecule (tetraeder of hydrogens)
+  TestMolecule = World::getInstance().createMolecule();
+  CPPUNIT_ASSERT(TestMolecule != NULL && "could not create molecule");
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1.5, 0., 1.5 ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1.5, 1.5 ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1.5, 1.5, 0. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0.5, 0.5, 0.5 ));
+  TestMolecule->AddAtom(Walker);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 5 );
+
+  // create stream with table
+  std::stringstream test;
+  test << ".\tH\tHe\tLi\tBe\tB\tC\n";
+  test << "H\t1.\t1.\t1.\t1.\t1.\t1.2\n";
+  test << "He\t1.\t1.\t1.\t1.\t1.\t1.\n";
+  test << "Li\t1.\t1.\t1.\t1.\t1.\t1.\n";
+  test << "Be\t1.\t1.\t1.\t1.\t1.\t1.\n";
+  test << "B\t1.\t1.\t1.\t1.\t1.\t1.\n";
+  test << "C\t1.2\t1.\t1.\t1.\t1.\t1.5\n";
+  BG = new BondGraph(true);
+  CPPUNIT_ASSERT(BG != NULL && "could not create BondGraph");
+
+  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(test) );
+  CPPUNIT_ASSERT_EQUAL( 1., BG->GetBondLength(0,0) );
+  CPPUNIT_ASSERT_EQUAL( 1.2, BG->GetBondLength(0,5) );
+  CPPUNIT_ASSERT_EQUAL( 1.5, BG->GetBondLength(5,5) );
+
+  molecule::atomVector Set = TestMolecule->getAtomSet();
+  BG->CreateAdjacency(Set);
+};
+
+
+void AnalysisBondsTest::tearDown()
+{
+  // remove the file
+  delete(BG);
+
+  // remove molecule
+  World::getInstance().destroyMolecule(TestMolecule);
+  // note that all the atoms are cleaned by TestMolecule
+  World::purgeInstance();
+};
+
+/** UnitTest for GetMaxMinMeanBondCount().
+ */
+void AnalysisBondsTest::GetMaxMinMeanBondCountTest()
+{
+  double Min = 20.; // check that initialization resets these arbitrary values
+  double Mean = 200.;
+  double Max = 1e-6;
+  GetMaxMinMeanBondCount(TestMolecule, Min, Mean, Max);
+  CPPUNIT_ASSERT_EQUAL( 1., Min );
+  CPPUNIT_ASSERT_EQUAL( 1.6, Mean );
+  CPPUNIT_ASSERT_EQUAL( 4., Max );
+
+};
+
+/** UnitTest for MinMaxBondDistanceBetweenElements().
+ */
+void AnalysisBondsTest::MinMeanMaxBondDistanceBetweenElementsTest()
+{
+  double Min = 20.; // check that initialization resets these arbitrary values
+  double Mean = 2e+6;
+  double Max = 1e-6;
+  double Min2 = 20.;
+  double Mean2 = 2e+6;
+  double Max2 = 1e-6;
+  const double maxbondlength = sqrt(1.*1. + 1.*1. + .5*.5);
+  const double minbondlength = sqrt(.5*.5 + .5*.5 + .5*.5);
+  const double meanbondlength = (minbondlength+3.*maxbondlength)/4.;
+  // check bond lengths C-H
+  MinMeanMaxBondDistanceBetweenElements(TestMolecule, hydrogen, carbon, Min, Mean, Max);
+  CPPUNIT_ASSERT_EQUAL( minbondlength , Min );
+  CPPUNIT_ASSERT_EQUAL( meanbondlength , Mean );
+  CPPUNIT_ASSERT_EQUAL( maxbondlength , Max );
+
+  // check that elements are symmetric, i.e. C-H == H-C
+  MinMeanMaxBondDistanceBetweenElements(TestMolecule, carbon, hydrogen, Min2, Mean2, Max2);
+  CPPUNIT_ASSERT_EQUAL( Min , Min2 );
+  CPPUNIT_ASSERT_EQUAL( Mean , Mean2 );
+  CPPUNIT_ASSERT_EQUAL( Max , Max2 );
+
+  // check no bond case (no bonds H-H in system!)
+  MinMeanMaxBondDistanceBetweenElements(TestMolecule, hydrogen, hydrogen, Min, Mean, Max);
+  CPPUNIT_ASSERT_EQUAL( 0. , Min );
+  CPPUNIT_ASSERT_EQUAL( 0. , Mean );
+  CPPUNIT_ASSERT_EQUAL( 0. , Max );
+};
Index: src/Analysis/unittests/AnalysisBondsUnitTest.hpp
===================================================================
--- src/Analysis/unittests/AnalysisBondsUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisBondsUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,49 @@
+/*
+ * AnalysisBondsUnitTest.hpp
+ *
+ *  Created on: Nov 7, 2009
+ *      Author: heber
+ */
+
+#ifndef ANALYSISBONDSUNITTEST_HPP_
+#define ANALYSISBONDSUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class BondGraph;
+class element;
+class molecule;
+class periodentafel;
+
+/********************************************** Test classes **************************************/
+
+class AnalysisBondsTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( AnalysisBondsTest) ;
+    CPPUNIT_TEST ( GetMaxMinMeanBondCountTest );
+    CPPUNIT_TEST ( MinMeanMaxBondDistanceBetweenElementsTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void GetMaxMinMeanBondCountTest();
+      void MinMeanMaxBondDistanceBetweenElementsTest();
+
+private:
+
+      molecule *TestMolecule;
+      const element *hydrogen;
+      const element *carbon;
+
+      BondGraph *BG;
+      string *filename;
+};
+
+#endif /* ANALYSISBONDSUNITTEST_HPP_ */
Index: src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,148 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * AnalysisCorrelationToPointUnitTest.cpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
+
+#include "Analysis/analysis_correlation.hpp"
+
+#include "Descriptors/MoleculeDescriptor.hpp"
+
+#include "atom.hpp"
+#include "element.hpp"
+#include "molecule.hpp"
+#include "linkedcell.hpp"
+#include "periodentafel.hpp"
+#include "World.hpp"
+
+#include "AnalysisCorrelationToPointUnitTest.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( AnalysisCorrelationToPointUnitTest );
+
+void AnalysisCorrelationToPointUnitTest::setUp()
+{
+  atom *Walker = NULL;
+
+  // init private all pointers to zero
+  TestMolecule = NULL;
+  pointmap = NULL;
+  binmap = NULL;
+  point = NULL;
+
+  // construct element list
+  std::vector<const element *> elements;
+  hydrogen = World::getInstance().getPeriode()->FindElement(1);
+  CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
+  elements.push_back(hydrogen);
+  // construct molecule (tetraeder of hydrogens)
+  TestMolecule = World::getInstance().createMolecule();
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
+  TestMolecule->AddAtom(Walker);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
+
+  TestMolecule->ActiveFlag = true;
+
+  // init point
+  point = new Vector(1.,1.,1.);
+
+  // init maps
+  World::getInstance().selectAllMolecules(AllMolecules());
+  allMolecules = World::getInstance().getSelectedMolecules();
+  CPPUNIT_ASSERT_EQUAL( (size_t) 1, allMolecules.size());
+  pointmap = CorrelationToPoint( allMolecules, elements, (const Vector *)point );
+  binmap = NULL;
+
+};
+
+
+void AnalysisCorrelationToPointUnitTest::tearDown()
+{
+  if (pointmap != NULL)
+    delete(pointmap);
+  if (binmap != NULL)
+    delete(binmap);
+
+  delete(point);
+  World::purgeInstance();
+  logger::purgeInstance();
+};
+
+void AnalysisCorrelationToPointUnitTest::CorrelationToPointTest()
+{
+  // do the pair correlation
+  CPPUNIT_ASSERT( pointmap != NULL );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, pointmap->size() );
+
+};
+
+void AnalysisCorrelationToPointUnitTest::CorrelationToPointBinNoRangeTest()
+{
+  BinPairMap::iterator tester;
+  // put pair correlation into bins and check with no range
+  binmap = BinData( pointmap, 0.5, 0., 0. );
+  OutputCorrelationMap<BinPairMap> ( (ofstream *)&cout, binmap, OutputCorrelation_Header, OutputCorrelation_Value );
+  CPPUNIT_ASSERT_EQUAL( (size_t)3, binmap->size() );
+  tester = binmap->begin();
+  CPPUNIT_ASSERT_EQUAL( 1., tester->first );
+  CPPUNIT_ASSERT_EQUAL( 3, tester->second );
+
+};
+
+void AnalysisCorrelationToPointUnitTest::CorrelationToPointBinRangeTest()
+{
+  BinPairMap::iterator tester;
+  // ... and check with [0., 2.] range
+  binmap = BinData( pointmap, 0.5, 0., 2. );
+  OutputCorrelationMap<BinPairMap> ( (ofstream *)&cout, binmap, OutputCorrelation_Header, OutputCorrelation_Value );
+  CPPUNIT_ASSERT_EQUAL( (size_t)5, binmap->size() );
+  tester = binmap->begin();
+  CPPUNIT_ASSERT_EQUAL( 0., tester->first );
+  tester = binmap->find(1.);
+  CPPUNIT_ASSERT_EQUAL( 1., tester->first );
+  CPPUNIT_ASSERT_EQUAL( 3, tester->second );
+
+};
Index: src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,52 @@
+/*
+ * AnalysisCorrelationToPointUnitTest.hpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+#ifndef ANALYSISCORRELATIONTOPOINTUNITTEST_HPP_
+#define ANALYSISCORRELATIONTOPOINTUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class element;
+class molecule;
+class MoleculeListClass;
+class Vector;
+
+/********************************************** Test classes **************************************/
+
+class AnalysisCorrelationToPointUnitTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( AnalysisCorrelationToPointUnitTest ) ;
+    CPPUNIT_TEST ( CorrelationToPointTest );
+    CPPUNIT_TEST ( CorrelationToPointBinNoRangeTest );
+    CPPUNIT_TEST ( CorrelationToPointBinRangeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void CorrelationToPointTest();
+      void CorrelationToPointBinNoRangeTest();
+      void CorrelationToPointBinRangeTest();
+
+private:
+
+      std::vector<molecule *> allMolecules;
+      molecule *TestMolecule;
+      const element *hydrogen;
+
+      CorrelationToPointMap *pointmap;
+      BinPairMap *binmap;
+      Vector *point;
+};
+
+#endif /* ANALYSISCORRELATIONTOPOINTUNITTEST_HPP_ */
Index: src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,250 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * AnalysisCorrelationToSurfaceUnitTest.cpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
+
+#include "Analysis/analysis_correlation.hpp"
+#include "atom.hpp"
+#include "Tesselation/boundary.hpp"
+#include "CodePatterns/Assert.hpp"
+#include "Descriptors/MoleculeDescriptor.hpp"
+#include "element.hpp"
+#include "molecule.hpp"
+#include "linkedcell.hpp"
+#include "periodentafel.hpp"
+#include "PointCloudAdaptor.hpp"
+#include "Tesselation/tesselation.hpp"
+#include "World.hpp"
+
+#include "AnalysisCorrelationToSurfaceUnitTest.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( AnalysisCorrelationToSurfaceUnitTest );
+
+void AnalysisCorrelationToSurfaceUnitTest::setUp()
+{
+  ASSERT_DO(Assert::Throw);
+
+  setVerbosity(5);
+
+  atom *Walker = NULL;
+
+  // init private all pointers to zero
+  TestSurfaceMolecule = NULL;
+  surfacemap = NULL;
+  binmap = NULL;
+  Surface = NULL;
+  LC = NULL;
+
+  // prepare element list
+  hydrogen = World::getInstance().getPeriode()->FindElement(1);
+  CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
+  elements.clear();
+
+  // construct molecule (tetraeder of hydrogens) base
+  TestSurfaceMolecule = World::getInstance().createMolecule();
+
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule->getAtomCount(), 4 );
+
+  TestSurfaceMolecule->ActiveFlag = true;
+
+  // init tesselation and linked cell
+  Surface = new Tesselation;
+  PointCloudAdaptor<molecule> cloud(TestSurfaceMolecule, TestSurfaceMolecule->name);
+  LC = new LinkedCell(cloud, 5.);
+  FindNonConvexBorder(TestSurfaceMolecule, Surface, (const LinkedCell *&)LC, 2.5, NULL);
+
+  // add outer atoms
+  carbon = World::getInstance().getPeriode()->FindElement(6);
+  TestSurfaceMolecule = World::getInstance().createMolecule();
+  Walker = World::getInstance().createAtom();
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(4., 0., 4. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  Walker = World::getInstance().createAtom();
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0., 4., 4. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  Walker = World::getInstance().createAtom();
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(4., 4., 0. ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  // add inner atoms
+  Walker = World::getInstance().createAtom();
+  Walker->setType(carbon);
+  Walker->setPosition(Vector(0.5, 0.5, 0.5 ));
+  TestSurfaceMolecule->AddAtom(Walker);
+
+  World::getInstance().selectAllMolecules(AllMolecules());
+  allMolecules = World::getInstance().getSelectedMolecules();
+  CPPUNIT_ASSERT_EQUAL( (size_t) 2, allMolecules.size());
+
+  // init maps
+  surfacemap = NULL;
+  binmap = NULL;
+
+};
+
+
+void AnalysisCorrelationToSurfaceUnitTest::tearDown()
+{
+  if (surfacemap != NULL)
+    delete(surfacemap);
+  if (binmap != NULL)
+    delete(binmap);
+
+  delete(Surface);
+  // note that all the atoms are cleaned by TestMolecule
+  delete(LC);
+  World::purgeInstance();
+  logger::purgeInstance();
+};
+
+
+/** Checks whether setup() does the right thing.
+ */
+void AnalysisCorrelationToSurfaceUnitTest::SurfaceTest()
+{
+  CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule->getAtomCount() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)2, allMolecules.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)6, Surface->LinesOnBoundary.size() );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->TrianglesOnBoundary.size() );
+};
+
+void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceTest()
+{
+  // do the pair correlation
+  elements.push_back(hydrogen);
+  surfacemap = CorrelationToSurface( allMolecules, elements, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
+  CPPUNIT_ASSERT( surfacemap != NULL );
+  CPPUNIT_ASSERT_EQUAL( (size_t)4, surfacemap->size() );
+};
+
+void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceHydrogenBinNoRangeTest()
+{
+  BinPairMap::iterator tester;
+  elements.push_back(hydrogen);
+  surfacemap = CorrelationToSurface( allMolecules, elements, Surface, LC );
+  // put pair correlation into bins and check with no range
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
+  binmap = BinData( surfacemap, 0.5, 0., 0. );
+  CPPUNIT_ASSERT_EQUAL( (size_t)1, binmap->size() );
+  OutputCorrelationMap<BinPairMap> ( (ofstream *)&cout, binmap, OutputCorrelation_Header, OutputCorrelation_Value );
+  tester = binmap->begin();
+  CPPUNIT_ASSERT_EQUAL( 0., tester->first );
+  CPPUNIT_ASSERT_EQUAL( 4, tester->second );
+
+};
+
+void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceHydrogenBinRangeTest()
+{
+  BinPairMap::iterator tester;
+  elements.push_back(hydrogen);
+  surfacemap = CorrelationToSurface( allMolecules, elements, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
+  // ... and check with [0., 2.] range
+  binmap = BinData( surfacemap, 0.5, 0., 2. );
+  CPPUNIT_ASSERT_EQUAL( (size_t)5, binmap->size() );
+//  OutputCorrelation ( (ofstream *)&cout, binmap );
+  tester = binmap->begin();
+  CPPUNIT_ASSERT_EQUAL( 0., tester->first );
+  CPPUNIT_ASSERT_EQUAL( 4, tester->second );
+  tester = binmap->find(1.);
+  CPPUNIT_ASSERT_EQUAL( 1., tester->first );
+  CPPUNIT_ASSERT_EQUAL( 0, tester->second );
+
+};
+
+void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceCarbonBinNoRangeTest()
+{
+  BinPairMap::iterator tester;
+  elements.push_back(carbon);
+  surfacemap = CorrelationToSurface( allMolecules, elements, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
+  // put pair correlation into bins and check with no range
+  binmap = BinData( surfacemap, 0.5, 0., 0. );
+  //OutputCorrelation ( (ofstream *)&cout, binmap );
+  CPPUNIT_ASSERT_EQUAL( (size_t)9, binmap->size() );
+  // inside point is first and must have negative value
+  tester = binmap->lower_bound(4.25-0.5); // start depends on the min value and
+  CPPUNIT_ASSERT( tester != binmap->end() );
+  CPPUNIT_ASSERT_EQUAL( 3, tester->second );
+  // inner point
+  tester = binmap->lower_bound(0.);
+  CPPUNIT_ASSERT( tester != binmap->end() );
+  CPPUNIT_ASSERT_EQUAL( 1, tester->second );
+};
+
+void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceCarbonBinRangeTest()
+{
+  BinPairMap::iterator tester;
+  elements.push_back(carbon);
+  surfacemap = CorrelationToSurface( allMolecules, elements, Surface, LC );
+//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
+  // ... and check with [0., 2.] range
+  binmap = BinData( surfacemap, 0.5, -2., 4. );
+  //OutputCorrelation ( (ofstream *)&cout, binmap );
+  CPPUNIT_ASSERT_EQUAL( (size_t)13, binmap->size() );
+  // three outside points
+  tester = binmap->lower_bound(4.25-0.5);
+  CPPUNIT_ASSERT( tester != binmap->end() );
+  CPPUNIT_ASSERT_EQUAL( 3, tester->second );
+  // inner point
+  tester = binmap->lower_bound(0.);
+  CPPUNIT_ASSERT( tester != binmap->end() );
+  CPPUNIT_ASSERT_EQUAL( 1, tester->second );
+};
Index: src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp
===================================================================
--- src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,63 @@
+/*
+ * AnalysisCorrelationToSurfaceUnitTest.hpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+#ifndef ANALYSISCORRELATIONTOSURFACEUNITTEST_HPP_
+#define ANALYSISCORRELATIONTOSURFACEUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class element;
+class LinkedCell;
+class molecule;
+class MoleculeListClass;
+class periodentafel;
+class Tesselation;
+
+/********************************************** Test classes **************************************/
+
+class AnalysisCorrelationToSurfaceUnitTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( AnalysisCorrelationToSurfaceUnitTest ) ;
+    CPPUNIT_TEST ( SurfaceTest );
+    CPPUNIT_TEST ( CorrelationToSurfaceTest );
+    CPPUNIT_TEST ( CorrelationToSurfaceHydrogenBinNoRangeTest );
+    CPPUNIT_TEST ( CorrelationToSurfaceHydrogenBinRangeTest );
+    CPPUNIT_TEST ( CorrelationToSurfaceCarbonBinNoRangeTest );
+    CPPUNIT_TEST ( CorrelationToSurfaceCarbonBinRangeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void SurfaceTest();
+      void CorrelationToSurfaceTest();
+      void CorrelationToSurfaceHydrogenBinNoRangeTest();
+      void CorrelationToSurfaceHydrogenBinRangeTest();
+      void CorrelationToSurfaceCarbonBinNoRangeTest();
+      void CorrelationToSurfaceCarbonBinRangeTest();
+
+private:
+
+      std::vector<molecule *> allMolecules;
+      molecule *TestSurfaceMolecule;
+      const element *hydrogen;
+      const element *carbon;
+      std::vector<const element *> elements;
+
+      CorrelationToSurfaceMap *surfacemap;
+      BinPairMap *binmap;
+      Tesselation *Surface;
+      LinkedCell *LC;
+};
+
+#endif /* ANALYSISCORRELATIONTOSURFACEUNITTEST_HPP_ */
Index: src/Analysis/unittests/AnalysisPairCorrelationUnitTest.cpp
===================================================================
--- src/Analysis/unittests/AnalysisPairCorrelationUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisPairCorrelationUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,145 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * AnalysisPairCorrelationUnitTest.cpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
+
+#include "Analysis/analysis_correlation.hpp"
+#include "Descriptors/MoleculeDescriptor.hpp"
+
+#include "World.hpp"
+#include "atom.hpp"
+#include "Tesselation/boundary.hpp"
+#include "element.hpp"
+#include "molecule.hpp"
+#include "linkedcell.hpp"
+#include "periodentafel.hpp"
+#include "Tesselation/tesselation.hpp"
+#include "World.hpp"
+
+#include "AnalysisPairCorrelationUnitTest.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( AnalysisPairCorrelationUnitTest );
+
+
+void AnalysisPairCorrelationUnitTest::setUp()
+{
+  atom *Walker = NULL;
+
+  // init private all pointers to zero
+  TestMolecule = NULL;
+  correlationmap = NULL;
+  binmap = NULL;
+
+  // construct element list
+  std::vector<const element *> elements;
+  hydrogen = World::getInstance().getPeriode()->FindElement(1);
+  CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
+  elements.push_back(hydrogen);
+  elements.push_back(hydrogen);
+
+  // construct molecule (tetraeder of hydrogens)
+  TestMolecule = World::getInstance().createMolecule();
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 0., 1. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 1., 1. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(1., 1., 0. ));
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0., 0., 0. ));
+  TestMolecule->AddAtom(Walker);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
+
+  // init maps
+  World::getInstance().selectAllMolecules(AllMolecules());
+  allMolecules = World::getInstance().getSelectedMolecules();
+  CPPUNIT_ASSERT_EQUAL( (size_t) 1, allMolecules.size());
+  correlationmap = PairCorrelation( allMolecules, elements);
+  binmap = NULL;
+
+};
+
+
+void AnalysisPairCorrelationUnitTest::tearDown()
+{
+  if (correlationmap != NULL)
+    delete(correlationmap);
+  if (binmap != NULL)
+    delete(binmap);
+
+  // note that all the atoms are cleaned by TestMolecule
+  World::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
+};
+
+void AnalysisPairCorrelationUnitTest::PairCorrelationTest()
+{
+  // do the pair correlation
+  CPPUNIT_ASSERT( correlationmap != NULL );
+  CPPUNIT_ASSERT_EQUAL( (size_t)6, correlationmap->size() );
+};
+
+void AnalysisPairCorrelationUnitTest::PairCorrelationBinNoRangeTest()
+{
+  BinPairMap::iterator tester;
+  // put pair correlation into bins and check with no range
+  binmap = BinData( correlationmap, 0.5, 0., 0. );
+  CPPUNIT_ASSERT_EQUAL( (size_t)1, binmap->size() );
+  //OutputCorrelation ( binmap );
+  tester = binmap->begin();
+  CPPUNIT_ASSERT_EQUAL( sqrt(2.), tester->first );
+  CPPUNIT_ASSERT_EQUAL( 6, tester->second );
+};
+
+void AnalysisPairCorrelationUnitTest::PairCorrelationBinRangeTest()
+{
+  BinPairMap::iterator tester;
+  // ... and check with [0., 2.] range
+  binmap = BinData( correlationmap, 0.5, 0., 2. );
+  CPPUNIT_ASSERT_EQUAL( (size_t)5, binmap->size() );
+  //OutputCorrelation ( binmap );
+  tester = binmap->begin();
+  CPPUNIT_ASSERT_EQUAL( 0., tester->first );
+  tester = binmap->find(1.);
+  CPPUNIT_ASSERT_EQUAL( 1., tester->first );
+  CPPUNIT_ASSERT_EQUAL( 6, tester->second );
+};
Index: src/Analysis/unittests/AnalysisPairCorrelationUnitTest.hpp
===================================================================
--- src/Analysis/unittests/AnalysisPairCorrelationUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/AnalysisPairCorrelationUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,51 @@
+/*
+ * AnalysisPairCorrelationUnitTest.hpp
+ *
+ *  Created on: Oct 13, 2009
+ *      Author: heber
+ */
+
+#ifndef ANALYSISPAIRCORRELATIONUNITTEST_HPP_
+#define ANALYSISPAIRCORRELATIONUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class element;
+class molecule;
+class MoleculeListClass;
+class Vector;
+
+/********************************************** Test classes **************************************/
+
+class AnalysisPairCorrelationUnitTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( AnalysisPairCorrelationUnitTest ) ;
+    CPPUNIT_TEST ( PairCorrelationTest );
+    CPPUNIT_TEST ( PairCorrelationBinNoRangeTest );
+    CPPUNIT_TEST ( PairCorrelationBinRangeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void PairCorrelationTest();
+      void PairCorrelationBinNoRangeTest();
+      void PairCorrelationBinRangeTest();
+
+private:
+
+      std::vector<molecule *> allMolecules;
+      molecule *TestMolecule;
+      const element *hydrogen;
+
+      PairCorrelationMap *correlationmap;
+      BinPairMap *binmap;
+};
+
+#endif /* ANALYSISPAIRCORRELATIONUNITTEST_HPP_ */
Index: src/Analysis/unittests/CountBondsUnitTest.cpp
===================================================================
--- src/Analysis/unittests/CountBondsUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/CountBondsUnitTest.cpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,243 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
+/*
+ * CountBondsUnitTest.cpp
+ *
+ *  Created on: Mar 30, 2010
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <iostream>
+#include <stdio.h>
+#include <cstring>
+
+#include "CodePatterns/Assert.hpp"
+
+#include "Analysis/analysis_bonds.hpp"
+#include "atom.hpp"
+#include "Bond/bond.hpp"
+#include "element.hpp"
+#include "Graph/BondGraph.hpp"
+#include "molecule.hpp"
+#include "periodentafel.hpp"
+#include "World.hpp"
+
+#include "CountBondsUnitTest.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( CountBondsTest );
+
+
+void CountBondsTest::setUp()
+{
+  atom *Walker = NULL;
+
+  // construct element
+  hydrogen = World::getInstance().getPeriode()->FindElement(1);
+  oxygen = World::getInstance().getPeriode()->FindElement(8);
+  CPPUNIT_ASSERT(hydrogen != NULL && "could not find element hydrogen");
+  CPPUNIT_ASSERT(oxygen != NULL && "could not find element oxygen");
+
+  // construct molecule (water molecule)
+  TestMolecule1 = World::getInstance().createMolecule();
+  CPPUNIT_ASSERT(TestMolecule1 != NULL && "could not create first molecule");
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(-0.2418, 0.9350, 0. ));
+  TestMolecule1->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0.9658, 0., 0. ));
+  TestMolecule1->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(oxygen);
+  Walker->setPosition(Vector(0., 0., 0. ));
+  TestMolecule1->AddAtom(Walker);
+
+  TestMolecule2 = World::getInstance().createMolecule();
+  CPPUNIT_ASSERT(TestMolecule2 != NULL && "could not create second molecule");
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(-0.2418, 0.9350, 0. ));
+  TestMolecule2->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(hydrogen);
+  Walker->setPosition(Vector(0.9658, 0., 0. ));
+  TestMolecule2->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
+  Walker->setType(oxygen);
+  Walker->setPosition(Vector(0., 0., 0. ));
+  TestMolecule2->AddAtom(Walker);
+
+  molecules = World::getInstance().getMolecules();
+  CPPUNIT_ASSERT(molecules != NULL && "could not obtain list of molecules");
+  molecules->insert(TestMolecule1);
+  molecules->insert(TestMolecule2);
+
+  // check that TestMolecule was correctly constructed
+  CPPUNIT_ASSERT_EQUAL( TestMolecule1->getAtomCount(), 3 );
+  CPPUNIT_ASSERT_EQUAL( TestMolecule2->getAtomCount(), 3 );
+
+  // create a small file with table
+  BG = new BondGraph(true);
+  CPPUNIT_ASSERT(BG != NULL && "could not create BondGraph");
+
+  // construct bond graphs
+  molecule::atomVector Set1 = TestMolecule1->getAtomSet();
+  BG->CreateAdjacency(Set1);
+  CPPUNIT_ASSERT( TestMolecule1->getBondCount() != 0);
+  molecule::atomVector Set2 = TestMolecule2->getAtomSet();
+  BG->CreateAdjacency(Set2);
+  CPPUNIT_ASSERT( TestMolecule2->getBondCount() != 0);
+//  TestMolecule1->Output((ofstream *)&cout);
+//  TestMolecule1->OutputBondsList(std::cout);
+};
+
+
+void CountBondsTest::tearDown()
+{
+  // remove the file
+  delete(BG);
+
+  World::purgeInstance();
+};
+
+/** UnitTest for CountBondsTest::BondsOfTwoTest().
+ */
+void CountBondsTest::BondsOfTwoTest()
+{
+  CPPUNIT_ASSERT_EQUAL( 4 , CountBondsOfTwo(molecules, hydrogen, oxygen) );
+  CPPUNIT_ASSERT_EQUAL( 0 , CountBondsOfTwo(molecules, hydrogen, hydrogen) );
+  CPPUNIT_ASSERT_EQUAL( 0 , CountBondsOfTwo(molecules, oxygen, oxygen) );
+};
+
+/** UnitTest for CountBondsTest::BondsOfThreeTest().
+ */
+void CountBondsTest::BondsOfThreeTest()
+{
+  CPPUNIT_ASSERT_EQUAL( 2 , CountBondsOfThree(molecules, hydrogen, oxygen, hydrogen) );
+  CPPUNIT_ASSERT_EQUAL( 0 , CountBondsOfThree(molecules, oxygen, hydrogen, oxygen) );
+};
+
+void OutputTestMolecule(molecule *mol, const char *name)
+{
+  ofstream output(name);
+  mol->OutputXYZ(&output);
+  output.close();
+}
+
+/** UnitTest for CountBondsTest::HydrogenBridgeBondsTest().
+ */
+void CountBondsTest::HydrogenBridgeBondsTest()
+{
+  double *mirror = new double[3];
+  CPPUNIT_ASSERT(mirror != NULL && "could not create array of doubles");
+  for (int i=0;i<3;i++)
+    mirror[i] = -1.;
+  Vector Translator;
+
+  //OutputTestMolecule(TestMolecule1, "testmolecule1.xyz");
+
+  cout << "Case 1: offset of (3,0,0), hence angles are (104.5, 0, 75.5, 180) < 30." << endl;
+  Translator  = Vector(3,0,0);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, oxygen, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-1.xyz");
+  Translator = Vector(-3,0,0);
+  TestMolecule2->Translate(&Translator);
+
+  cout << "Case 2: offset of (0,3,0), hence angle are (14.5, 165.5, 90) < 30 (only three, because other 90 is missing due to first H01 only fulfilling H-bond criteria)." << endl;
+  Translator = Vector(0,3,0);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-2.xyz");
+  Translator = Vector(0,-3,0);
+  TestMolecule2->Translate(&Translator);
+
+  cout << "Case 3: offset of (0,-3,0) and mirror, hence angle are (165.5, 90, 165.5, 90) > 30." << endl;
+  Translator = Vector(0,-3,0);
+  TestMolecule2->Scale((const double ** const)&mirror);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-3.xyz");
+  Translator = Vector(0,3,0);
+  TestMolecule2->Translate(&Translator);
+  TestMolecule2->Scale((const double ** const)&mirror);
+
+  cout << "Case 4: offset of (2,1,0), hence angle are (78, 26.6, 102, 153.4) < 30." << endl;
+  Translator = Vector(2,1,0);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-4.xyz");
+  Translator = Vector(-2,-1,0);
+  TestMolecule2->Translate(&Translator);
+
+  cout << "Case 5: offset of (0,0,3), hence angle are (90, 90, 90, 90) > 30." << endl;
+  Translator = Vector(0,0,3);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-5.xyz");
+  Translator = Vector(0,0,-3);
+  TestMolecule2->Translate(&Translator);
+
+  cout << "Case 6: offset of (-3,0,0) and mirror, hence angle are (75.5, 180, 104.5, 180) > 30." << endl;
+  Translator = Vector(-3,0,0);
+  TestMolecule2->Scale((const double ** const)&mirror);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-6.xyz");
+  Translator = Vector(3,0,0);
+  TestMolecule2->Translate(&Translator);
+  TestMolecule2->Scale((const double ** const)&mirror);
+
+  cout << "Case 7: offset of (3,0,0) and mirror, hence angles are (104.5, 0, 104.5, 0) < 30, but interfering hydrogens." << endl;
+  Translator = Vector(3,0,0);
+  TestMolecule2->Scale((const double ** const)&mirror);
+  TestMolecule2->Translate(&Translator);
+  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-7.xyz");
+  Translator = Vector(-3,0,0);
+  TestMolecule2->Translate(&Translator);
+  TestMolecule2->Scale((const double ** const)&mirror);
+
+  cout << "Case 8: offset of (0,3,0), hence angle are (14.5, 90, 14.5, 90) < 30, but interfering hydrogens." << endl;
+  Translator = Vector(0,3,0);
+  TestMolecule2->Scale((const double ** const)&mirror);
+  TestMolecule2->Translate(&Translator);
+  //OutputTestMolecule(TestMolecule2, "testmolecule2-8.xyz");
+  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
+  Translator = Vector(0,-3,0);
+  TestMolecule2->Translate(&Translator);
+  TestMolecule2->Scale((const double ** const)&mirror);
+
+  delete[](mirror);
+};
Index: src/Analysis/unittests/CountBondsUnitTest.hpp
===================================================================
--- src/Analysis/unittests/CountBondsUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/CountBondsUnitTest.hpp	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,55 @@
+/*
+ * CountBondsUnitTest.hpp
+ *
+ *  Created on: Mar 30, 2010
+ *      Author: heber
+ */
+
+#ifndef COUNTBONDSUNITTEST_HPP_
+#define COUNTBONDSUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <cppunit/extensions/HelperMacros.h>
+
+
+class BondGraph;
+class element;
+class molecule;
+class periodentafel;
+
+/********************************************** Test classes **************************************/
+
+class CountBondsTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( CountBondsTest) ;
+    CPPUNIT_TEST ( BondsOfTwoTest );
+    CPPUNIT_TEST ( BondsOfThreeTest );
+    CPPUNIT_TEST ( HydrogenBridgeBondsTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void BondsOfTwoTest();
+      void BondsOfThreeTest();
+      void HydrogenBridgeBondsTest();
+
+private:
+
+      MoleculeListClass *molecules;
+      molecule *TestMolecule1;
+      molecule *TestMolecule2;
+      const element *hydrogen;
+      const element *oxygen;
+
+      BondGraph *BG;
+      string *filename;
+};
+
+
+#endif /* COUNTBONDSUNITTEST_HPP_ */
Index: src/Analysis/unittests/Makefile.am
===================================================================
--- src/Analysis/unittests/Makefile.am	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
+++ src/Analysis/unittests/Makefile.am	(revision 9b5a2c623129c569571d4a77eabff8bb2bc84231)
@@ -0,0 +1,63 @@
+# PLEASE adhere to the alphabetical ordering in this Makefile!
+# Also indentation by a single tab
+
+ANALYSISTESTSSOURCES = \
+  ../Analysis/unittests/AnalysisBondsUnitTest.cpp \
+  ../Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp \
+  ../Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp  \
+  ../Analysis/unittests/AnalysisPairCorrelationUnitTest.cpp \
+  ../Analysis/unittests/CountBondsUnitTest.cpp
+
+ANALYSISTESTSHEADERS = \
+  ../Analysis/unittests/AnalysisBondsUnitTest.hpp \
+  ../Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp \
+  ../Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp  \
+  ../Analysis/unittests/AnalysisPairCorrelationUnitTest.hpp \
+  ../Analysis/unittests/CountBondsUnitTest.hpp
+
+ANALYSISTESTS = \
+  AnalysisBondsUnitTest \
+  AnalysisCorrelationToPointUnitTest \
+  AnalysisCorrelationToSurfaceUnitTest \
+  AnalysisPairCorrelationUnitTest \
+  CountBondsUnitTest
+
+TESTS += $(ANALYSISTESTS)
+check_PROGRAMS += $(ANALYSISTESTS) 
+noinst_PROGRAMS += $(ANALYSISTESTS)
+
+ANALYSISLIBS = \
+	../libMolecuilderAnalysis.la \
+	../libMolecuilderUI.la \
+	../libMolecuilder.la \
+	$(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
+	${CodePatterns_LIBS} \
+	$(BOOST_LIB)
+
+AnalysisBondsUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Analysis/unittests/AnalysisBondsUnitTest.cpp \
+	../Analysis/unittests/AnalysisBondsUnitTest.hpp
+AnalysisBondsUnitTest_LDADD = ${ANALYSISLIBS}
+
+AnalysisCorrelationToPointUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Analysis/unittests/AnalysisCorrelationToPointUnitTest.cpp \
+	../Analysis/unittests/AnalysisCorrelationToPointUnitTest.hpp
+AnalysisCorrelationToPointUnitTest_LDADD = ${ANALYSISLIBS}
+
+AnalysisCorrelationToSurfaceUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp \
+	../Analysis/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp
+AnalysisCorrelationToSurfaceUnitTest_LDADD = ${ANALYSISLIBS}
+
+AnalysisPairCorrelationUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Analysis/unittests/AnalysisPairCorrelationUnitTest.cpp \
+	../Analysis/unittests/AnalysisPairCorrelationUnitTest.hpp
+AnalysisPairCorrelationUnitTest_LDADD = ${ANALYSISLIBS}
+
+CountBondsUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../Analysis/unittests/CountBondsUnitTest.cpp \
+	../Analysis/unittests/CountBondsUnitTest.hpp
+CountBondsUnitTest_LDADD = ${ANALYSISLIBS}
+
+
+#AUTOMAKE_OPTIONS = parallel-tests
