Index: molecuilder/src/analysis_correlation.cpp
===================================================================
--- molecuilder/src/analysis_correlation.cpp	(revision 20895b5069ce1fb4aba27283ed8d9de1971e75fc)
+++ molecuilder/src/analysis_correlation.cpp	(revision 5f1d0216136640784700cd053eed75cf132babd9)
@@ -65,4 +65,81 @@
 };
 
+/** Calculates the pair correlation between given elements.
+ * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
+ * \param *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type1 first element or NULL (if any element)
+ * \param *type2 second element or NULL (if any element)
+ * \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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
+{
+  PairCorrelationMap *outmap = NULL;
+  double distance = 0.;
+  int n[NDIM];
+  Vector checkX;
+  Vector periodicX;
+  int Othern[NDIM];
+  Vector checkOtherX;
+  Vector periodicOtherX;
+
+  if (molecules->ListOfMolecules.empty()) {
+    cerr << Verbose(1) <<"No molecule given." << endl;
+    return outmap;
+  }
+  outmap = new PairCorrelationMap;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
+      cerr << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+        if ((type1 == NULL) || (Walker->type == type1)) {
+          periodicX.CopyVector(Walker->node);
+          periodicX.MatrixMultiplication(FullInverseMatrix);  // 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.Init(n[0], n[1], n[2]);
+                checkX.AddVector(&periodicX);
+                checkX.MatrixMultiplication(FullMatrix);
+                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
+                  if ((*MolOtherWalker)->ActiveFlag) {
+                    *out << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
+                    atom *OtherWalker = (*MolOtherWalker)->start;
+                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
+                      OtherWalker = OtherWalker->next;
+                      *out << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
+                      if (Walker->nr < OtherWalker->nr)
+                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
+                          periodicOtherX.CopyVector(OtherWalker->node);
+                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // 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.Init(Othern[0], Othern[1], Othern[2]);
+                                checkOtherX.AddVector(&periodicOtherX);
+                                checkOtherX.MatrixMultiplication(FullMatrix);
+                                distance = checkX.Distance(&checkOtherX);
+                                //*out << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
+                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
+                              }
+                        }
+                  }
+              }
+            }
+        }
+      }
+    }
+
+  return outmap;
+};
+
 /** Calculates the distance (pair) correlation between a given element and a point.
  * \param *out output stream for debugging
@@ -100,4 +177,55 @@
 };
 
+/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
+ * \param *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type element or NULL (if any element)
+ * \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(  ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
+{
+  CorrelationToPointMap *outmap = NULL;
+  double distance = 0.;
+  int n[NDIM];
+  Vector periodicX;
+  Vector checkX;
+
+  if (molecules->ListOfMolecules.empty()) {
+    *out << Verbose(1) <<"No molecule given." << endl;
+    return outmap;
+  }
+  outmap = new CorrelationToPointMap;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
+      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+        if ((type == NULL) || (Walker->type == type)) {
+          periodicX.CopyVector(Walker->node);
+          periodicX.MatrixMultiplication(FullInverseMatrix);  // 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.Init(n[0], n[1], n[2]);
+                checkX.AddVector(&periodicX);
+                checkX.MatrixMultiplication(FullMatrix);
+                distance = checkX.Distance(point);
+                *out << Verbose(4) << "Current distance is " << distance << "." << endl;
+                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
+              }
+        }
+      }
+    }
+
+  return outmap;
+};
+
 /** Calculates the distance (pair) correlation between a given element and a surface.
  * \param *out output stream for debugging
@@ -114,5 +242,49 @@
   class BoundaryTriangleSet *triangle = NULL;
   Vector centroid;
-  int ranges[NDIM] = {1, 1, 1};
+
+  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
+    *out << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
+    return outmap;
+  }
+  outmap = new CorrelationToSurfaceMap;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+        if ((type == NULL) || (Walker->type == type)) {
+          triangle = Surface->FindClosestTriangleToPoint(out, Walker->node, LC );
+          if (triangle != NULL) {
+            distance = DistanceToTrianglePlane(out, Walker->node, triangle);
+            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, 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 *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type element or NULL (if any element)
+ * \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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
+{
+  CorrelationToSurfaceMap *outmap = NULL;
+  double distance = 0;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector centroid;
   int n[NDIM];
   Vector periodicX;
Index: molecuilder/src/analysis_correlation.hpp
===================================================================
--- molecuilder/src/analysis_correlation.hpp	(revision 20895b5069ce1fb4aba27283ed8d9de1971e75fc)
+++ molecuilder/src/analysis_correlation.hpp	(revision 5f1d0216136640784700cd053eed75cf132babd9)
@@ -48,4 +48,7 @@
 CorrelationToPointMap *CorrelationToPoint( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point );
 CorrelationToSurfaceMap *CorrelationToSurface( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC );
+PairCorrelationMap *PeriodicPairCorrelation( ofstream * const out, MoleculeListClass * const &molecules, const element * const type1, const element * const  type2, const int ranges[NDIM] );
+CorrelationToPointMap *PeriodicCorrelationToPoint( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] );
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
 double GetBin ( const double value, const double BinWidth, const double BinStart );
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map );
Index: molecuilder/src/builder.cpp
===================================================================
--- molecuilder/src/builder.cpp	(revision 20895b5069ce1fb4aba27283ed8d9de1971e75fc)
+++ molecuilder/src/builder.cpp	(revision 5f1d0216136640784700cd053eed75cf132babd9)
@@ -1632,5 +1632,6 @@
                 element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
                 FindNonConvexBorder((ofstream *)&cout, Boundary, TesselStruct, LCList, radius, NULL);
-                CorrelationToSurfaceMap *surfacemap = CorrelationToSurface( (ofstream *)&cout, molecules, elemental, TesselStruct, LCList );
+                int ranges[NDIM] = {1,1,1};
+                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( (ofstream *)&cout, molecules, elemental, TesselStruct, LCList, ranges );
                 BinPairMap *binmap = BinData( (ofstream *)&cout, surfacemap, 0.5, 0., 0. );
                 OutputCorrelation ( &binoutput, binmap );
Index: molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
===================================================================
--- molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 20895b5069ce1fb4aba27283ed8d9de1971e75fc)
+++ molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 5f1d0216136640784700cd053eed75cf132babd9)
@@ -186,5 +186,5 @@
   CPPUNIT_ASSERT_EQUAL( 3, tester->second );
   // inner point
-  tester = binmap->lower_bound(0.);
+  tester = binmap->lower_bound(-0.5);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 1, tester->second );
@@ -200,9 +200,9 @@
   OutputCorrelation ( (ofstream *)&cout, binmap );
   // three outside points
-  tester = binmap->find(3.);
+  tester = binmap->lower_bound(3.);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 3, tester->second );
   // inner point
-  tester = binmap->find(0.);
+  tester = binmap->lower_bound(-0.5);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 1, tester->second );
Index: molecuilder/src/vector.cpp
===================================================================
--- molecuilder/src/vector.cpp	(revision 20895b5069ce1fb4aba27283ed8d9de1971e75fc)
+++ molecuilder/src/vector.cpp	(revision 5f1d0216136640784700cd053eed75cf132babd9)
@@ -276,5 +276,8 @@
   temp.SubtractVector(PlaneOffset);
   double sign = temp.ScalarProduct(PlaneNormal);
-  sign /= fabs(sign);
+  if (fabs(sign) > MYEPSILON)
+    sign /= fabs(sign);
+  else
+    sign = 0.;
 
   return (temp.Norm()*sign);
