Index: src/analysis_correlation.cpp
===================================================================
--- src/analysis_correlation.cpp	(revision 164a33ff15a1bf3ee4cf3384e5372f87f6e826e0)
+++ src/analysis_correlation.cpp	(revision 775d1332dac2d98e150640bab72e75a8a55d0a66)
@@ -372,5 +372,5 @@
   *file << "BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << runner->first << "\t" << runner->second << endl;
+    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
   }
 };
@@ -385,5 +385,5 @@
   *file << "BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
   }
 };
@@ -400,5 +400,5 @@
     *file << runner->first;
     for (int i=0;i<NDIM;i++)
-      *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
+      *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
     *file << endl;
   }
@@ -414,6 +414,6 @@
   *file << "BinStart\tTriangle" << endl;
   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
-  }
-};
-
+    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+  }
+};
+
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 164a33ff15a1bf3ee4cf3384e5372f87f6e826e0)
+++ src/boundary.cpp	(revision 775d1332dac2d98e150640bab72e75a8a55d0a66)
@@ -789,4 +789,5 @@
  * \param *filler molecule which the box is to be filled with
  * \param configuration contains box dimensions
+ * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
  * \param distance[NDIM] distance between filling molecules in each direction
  * \param boundary length of boundary zone between molecule and filling mollecules
@@ -797,5 +798,5 @@
  * \return *mol pointer to new molecule with filled atoms
  */
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
 {
 	Info FunctionInfo(__func__);
@@ -862,9 +863,9 @@
           } else {
             const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
-            FillIt = FillIt && (distance > boundary*boundary);
+            FillIt = FillIt && (distance > boundary*boundary) && ((MaxDistance < 0) || (MaxDistance*MaxDistance > distance));
             if (FillIt) {
               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
             } else {
-              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
+              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point, within boundary or outside of MaxDistance." << endl;
               break;
             }
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 164a33ff15a1bf3ee4cf3384e5372f87f6e826e0)
+++ src/boundary.hpp	(revision 775d1332dac2d98e150640bab72e75a8a55d0a66)
@@ -49,5 +49,5 @@
 
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
 void FindConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 164a33ff15a1bf3ee4cf3384e5372f87f6e826e0)
+++ src/builder.cpp	(revision 775d1332dac2d98e150640bab72e75a8a55d0a66)
@@ -1387,4 +1387,5 @@
   enum ConfigStatus configPresent = absent;
   clock_t start,end;
+  double MaxDistance = -1;
   int argptr;
   molecule *mol = NULL;
@@ -1412,5 +1413,5 @@
             Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
             Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
-            Log() << Verbose(0) << "\t-C <type> [parms] <output> <bin output>\tPair Correlation analysis." << endl;
+            Log() << Verbose(0) << "\t-C <type> [params] <output> <bin output> <BinStart> <BinEnd>\tPair Correlation analysis." << endl;
             Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
@@ -1419,4 +1420,5 @@
             Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
             Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
+            Log() << Verbose(0) << "\t-FF <MaxDistance> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
             Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
             Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
@@ -1680,5 +1682,5 @@
               if ((argptr >= argc)) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl;
+                eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl;
                 performCriticalExit();
               } else {
@@ -1686,5 +1688,5 @@
                   case 'E':
                     {
-                      if ((argptr+4 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) {
+                      if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) {
                         ExitFlag = 255;
                         eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl;
@@ -1693,10 +1695,12 @@
                         ofstream output(argv[argptr+3]);
                         ofstream binoutput(argv[argptr+4]);
-
-                        element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
-                        element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+1]));
+                        const double BinStart = atof(argv[argptr+5]);
+                        const double BinEnd = atof(argv[argptr+6]);
+
+                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
+                        element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2]));
                         PairCorrelationMap *correlationmap = PairCorrelation(molecules, elemental, elemental2);
                         //OutputCorrelationToSurface(&output, correlationmap);
-                        BinPairMap *binmap = BinData( correlationmap, 0.5, 0., 0. );
+                        BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
                         OutputCorrelation ( &binoutput, binmap );
                         output.close();
@@ -1704,5 +1708,5 @@
                         delete(binmap);
                         delete(correlationmap);
-                        argptr+=5;
+                        argptr+=7;
                       }
                     }
@@ -1711,5 +1715,5 @@
                   case 'P':
                     {
-                      if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) {
+                      if ((argptr+8 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+7])) || (!IsValidNumber(argv[argptr+8])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) {
                         ExitFlag = 255;
                         eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl;
@@ -1718,10 +1722,12 @@
                         ofstream output(argv[argptr+5]);
                         ofstream binoutput(argv[argptr+6]);
-
-                        element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
+                        const double BinStart = atof(argv[argptr+7]);
+                        const double BinEnd = atof(argv[argptr+8]);
+
+                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
                         Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3]));
                         CorrelationToPointMap *correlationmap = CorrelationToPoint(molecules, elemental, Point);
                         //OutputCorrelationToSurface(&output, correlationmap);
-                        BinPairMap *binmap = BinData( correlationmap, 0.5, 0., 0. );
+                        BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
                         OutputCorrelation ( &binoutput, binmap );
                         output.close();
@@ -1730,5 +1736,5 @@
                         delete(binmap);
                         delete(correlationmap);
-                        argptr+=7;
+                        argptr+=9;
                       }
                     }
@@ -1737,5 +1743,5 @@
                   case 'S':
                     {
-                      if ((argptr+3 >= argc) || (!IsValidNumber(argv[argptr+1])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) {
+                      if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) {
                         ExitFlag = 255;
                         eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C S <Z> <output> <bin output>" << endl;
@@ -1745,5 +1751,6 @@
                         ofstream binoutput(argv[argptr+3]);
                         const double radius = 5.;
-
+                        const double BinStart = atof(argv[argptr+4]);
+                        const double BinEnd = atof(argv[argptr+5]);
 
                         // get the boundary
@@ -1766,10 +1773,10 @@
                         }
                         LCList = new LinkedCell(Boundary, 2.*radius);
-                        element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
+                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
                         FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
                         int ranges[NDIM] = {1,1,1};
                         CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
-                        //OutputCorrelationToSurface(&output, surfacemap);
-                        BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
+                        OutputCorrelationToSurface(&output, surfacemap);
+                        BinPairMap *binmap = BinData( surfacemap, 0.5, BinStart, BinEnd );
                         OutputCorrelation ( &binoutput, binmap );
                         output.close();
@@ -1782,5 +1789,5 @@
                         delete(binmap);
                         delete(surfacemap);
-                        argptr+=4;
+                        argptr+=6;
                       }
                     }
@@ -1811,4 +1818,10 @@
             case 'F':
               if (ExitFlag == 0) ExitFlag = 1;
+              MaxDistance = -1;
+              if (argv[argptr-1][2] == 'F') {
+                // fetch first argument as max distance to surface
+                MaxDistance = atof(argv[argptr++]);
+                Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl;
+              }
               if (argptr+6 >=argc) {
                 ExitFlag = 255;
@@ -1844,5 +1857,5 @@
                 for (int i=0;i<NDIM;i++)
                   distance[i] = atof(argv[argptr+i]);
-                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
+                Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
                 if (Filling != NULL) {
                   Filling->ActiveFlag = false;
