Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 110cebe3fea3969bcc2d7a77d60d8d02c0c7734f)
+++ src/boundary.cpp	(revision 6c58125d49933c6c85972c57b3af4f3c937fa82f)
@@ -169,4 +169,9 @@
 // ========================================== F U N C T I O N S =================================
 
+/** Finds the endpoint two lines are sharing.
+ * \param *line1 first line
+ * \param *line2 second line
+ * \return point which is shared or NULL if none
+ */
 class BoundaryPointSet * GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
 {
@@ -437,12 +442,14 @@
  * \param *out output stream for debugging
  * \param *configuration needed for path to store convex envelope file
+ * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
  * \param *mol molecule structure representing the cluster
  */
-double VolumeOfConvexEnvelope(ofstream *out, config *configuration, molecule *mol) 
+double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPtr, molecule *mol) 
 {
   bool IsAngstroem = configuration->GetIsAngstroem();
   atom *Walker = NULL;
   struct Tesselation *TesselStruct = new Tesselation;
-  Boundaries *BoundaryPoints = NULL;
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
   
   // 1. calculate center of gravity
@@ -459,5 +466,10 @@
   
   // 3. Find all points on the boundary
-  BoundaryPoints = GetBoundaryPoints(out, mol);
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
   
   // 3d. put into boundary set only those points appearing in each of the NDIM sets
@@ -544,5 +556,7 @@
 
   // free reference lists
-
+  if (BoundaryFreeFlag)
+    delete[](BoundaryPoints);
+  
   return volume;
 };
@@ -550,14 +564,23 @@
 
 /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
+ * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
  * \param *out output stream for debugging
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
- * \param clustervolume volume of the convex envelope cluster, \sa VolumeOfConvexEnvelope()
+ * \param repetition[] number of repeated cluster per axis
  * \param celldensity desired average density in final cell
- * \param GreatestDiameter[]  greatest diameter of the cluster, \sa GetDiametersOfCluster()
  */
-void CreateClustersinWater(ofstream *out, config *configuration, molecule *mol, double clustervolume, double celldensity, double GreatestDiameter[NDIM])
-{
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity)
+{
+  // some preparations beforehand
   bool IsAngstroem = configuration->GetIsAngstroem();
+  Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
+  double clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, IsAngstroem);
+  double coeffs[NDIM];
+  int TotalNoClusters = 1;
+  for (int i=0;i<NDIM;i++)
+    TotalNoClusters *= repetition[i];
+
   // sum up the atomic masses
   double totalmass = 0.;
@@ -575,17 +598,20 @@
   double cellvolume;
   if (IsAngstroem)
-    cellvolume = (4*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
   else
-    cellvolume = (4*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
   *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  double minimumvolume = 4*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
+  double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
   *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
   if (minimumvolume < cellvolume)
     cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
-  double  a = 4*(GreatestDiameter[0]+GreatestDiameter[1]+GreatestDiameter[2]);
-  double b = 4*(GreatestDiameter[0]*GreatestDiameter[1] + GreatestDiameter[0]*GreatestDiameter[2] + GreatestDiameter[1]*GreatestDiameter[2]);
-  double c = minimumvolume - cellvolume;
+  
+  coeffs[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
+  coeffs[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
+            + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2] 
+            + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
+  coeffs[2] = minimumvolume - cellvolume;
   double x0 = 0.,x1 = 0.,x2 = 0.;
-  if (gsl_poly_solve_cubic(a,b,c,&x0,&x1,&x2) == 1) // either 1 or 3 on return
+  if (gsl_poly_solve_cubic(coeffs[0],coeffs[1],coeffs[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
     *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
   else {
@@ -593,9 +619,11 @@
     x0 = x2;  // sorted in ascending order
   }
-  
-  a = 2 * (x0 + GreatestDiameter[0]);
-  b = 2 * (x0 + GreatestDiameter[1]);
-  c = (x0 + GreatestDiameter[2]);
-  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << a << " and " << b << " and " << c << " with total volume of " << a*b*c << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+
+  cellvolume = 1;
+  for(int i=0;i<NDIM;i++) {
+    coeffs[i] = repetition[0] * (x0 + GreatestDiameter[0]);
+    cellvolume *= coeffs[i];
+  }
+  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << coeffs[0] << " and " << coeffs[1] << " and " << coeffs[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
 };
 
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 110cebe3fea3969bcc2d7a77d60d8d02c0c7734f)
+++ src/boundary.hpp	(revision 6c58125d49933c6c85972c57b3af4f3c937fa82f)
@@ -122,6 +122,6 @@
 
 
-double VolumeOfConvexEnvelope(ofstream *out, config *configuration, molecule *mol);
-void CreateClustersinWater(ofstream *out, config *configuration, molecule *mol, double clustervolume, double celldensity, double GreatestDiameter[NDIM]);
+double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPoints, molecule *mol);
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity);
 
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 110cebe3fea3969bcc2d7a77d60d8d02c0c7734f)
+++ src/builder.cpp	(revision 6c58125d49933c6c85972c57b3af4f3c937fa82f)
@@ -567,5 +567,5 @@
     case 'e':
     	cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-    	VolumeOfConvexEnvelope((ofstream *)&cout, configuration, mol);
+    	VolumeOfConvexEnvelope((ofstream *)&cout, configuration, NULL, mol);
     	break;
   }
@@ -734,5 +734,4 @@
   string line;
   atom *first;
-  double tmp;
   int ExitFlag = 0;
   int j;
@@ -989,5 +988,5 @@
 //                tmp = 1.3;
 //              }
-              VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, mol);
+              VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
               break;
             default:   // no match? Step on
