Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 6c58125d49933c6c85972c57b3af4f3c937fa82f)
+++ src/boundary.cpp	(revision 318bfd37d3e5e3924d9bac94b3435c4b628d634e)
@@ -372,9 +372,20 @@
  * \param *out output stream for debugging
  * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
+ * \param *mol molecule structure representing the cluster
  * \param IsAngstroem whether we have angstroem or atomic units
  * \return NDIM array of the diameters
  */ 
-double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPoints, bool IsAngstroem)
-{
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
+{
+  // get points on boundary of NULL was given as parameter
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
+  
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
@@ -434,4 +445,8 @@
   *out << Verbose(0) << "RESULT: The biggest diameters are " << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and " << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "." << endl;
 
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[](BoundaryPoints);
+
   return GreatestDiameter;
 };
@@ -452,5 +467,10 @@
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
-  
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G,h;
+  vector x,y;
+  double a,b,c;
+
   // 1. calculate center of gravity
   *out << endl;
@@ -473,31 +493,9 @@
   }
   
-  // 3d. put into boundary set only those points appearing in each of the NDIM sets
-  int *AtomList = new int[mol->AtomCount];
-  for (int j=0; j<mol->AtomCount; j++) // reset list
-    AtomList[j] = 0;
-  for (int axis=0; axis<NDIM; axis++)  {  // fill list when it's on the boundary
+  // 4. fill the boundary point list
+  for (int axis=0;axis<NDIM;axis++)
     for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-      AtomList[runner->second.second->nr]++;
+      TesselStruct->AddPoint(runner->second.second);
     }
-  }
-  for (int axis=0; axis<NDIM; axis++)  {
-    for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-      if (AtomList[runner->second.second->nr] < 1) {
-        *out << Verbose(1) << "Throwing especially out " << *runner->second.second << " in axial projection of axis " << axis << "." << endl;
-        BoundaryPoints[axis].erase(runner);
-      }
-    }
-  }
-  delete[](AtomList);
-  
-  // 4a. fill the boundary point list
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    if (AtomList[Walker->nr] > 0) {
-      TesselStruct->AddPoint(Walker);
-    }
-  }
 
   *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
@@ -512,8 +510,8 @@
 //  *out << endl;
   
-  // 4b. guess starting triangle
+  // 5a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
   
-  // 5. go through all lines, that are not yet part of two triangles (only of one so far)
+  // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
   TesselStruct->TesselateOnBoundary(out, configuration, mol);
 
@@ -522,9 +520,4 @@
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
   *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
-  double volume = 0.;
-  double PyramidVolume = 0.;
-  double G,h;
-  vector x,y;
-  double a,b,c;
   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
     x.CopyVector(&runner->second->endpoints[0]->node->x);
@@ -568,15 +561,18 @@
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
- * \param repetition[] number of repeated cluster per axis
  * \param celldensity desired average density in final cell
  */
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity)
-{
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity)
+{
+  // transform to PAS
+  mol->PrincipalAxisSystem(out, true);
+  
   // 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];
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
+  vector BoxLengths;
+  int repetition[NDIM] = {1, 1, 1};
   int TotalNoClusters = 1;
   for (int i=0;i<NDIM;i++)
@@ -602,16 +598,17 @@
     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 = 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)
+  if (minimumvolume > cellvolume)
     cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
   
-  coeffs[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
-  coeffs[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
+  BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
+  BoxLengths.x[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;
+  BoxLengths.x[2] = minimumvolume - cellvolume;
   double x0 = 0.,x1 = 0.,x2 = 0.;
-  if (gsl_poly_solve_cubic(coeffs[0],coeffs[1],coeffs[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
+  if (gsl_poly_solve_cubic(BoxLengths.x[0],BoxLengths.x[1],BoxLengths.x[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
     *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
   else {
@@ -622,8 +619,15 @@
   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;
+    BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+    cellvolume *= BoxLengths.x[i];
+  }
+  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
+  // set new box dimensions
+  *out << Verbose(0) << "Translating to box with these boundaries." << endl;
+  mol->CenterInBox((ofstream *)&cout, &BoxLengths);
+  // update Box of atoms by boundary
+  mol->SetBoxDimension(&BoxLengths);
+  
 };
 
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 6c58125d49933c6c85972c57b3af4f3c937fa82f)
+++ src/boundary.hpp	(revision 318bfd37d3e5e3924d9bac94b3435c4b628d634e)
@@ -123,5 +123,6 @@
 
 double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPoints, molecule *mol);
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity);
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity);
 
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 6c58125d49933c6c85972c57b3af4f3c937fa82f)
+++ src/builder.cpp	(revision 318bfd37d3e5e3924d9bac94b3435c4b628d634e)
@@ -758,4 +758,5 @@
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
+            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
             cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
@@ -768,4 +769,5 @@
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
+            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 
             cout << "\t-v/-V\t\tGives version information." << endl;
             cout << "Note: config files must not begin with '-' !" << endl;
@@ -983,10 +985,76 @@
               ExitFlag = 1;
               cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-//              tmp = atof(argv[argptr++]);
-//              if (tmp < 1.0) {
-//                cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
-//                tmp = 1.3;
-//              }
               VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
+              break;
+            case 'u':
+              {
+                double tmp;
+                ExitFlag = 1;
+                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
+                tmp = atof(argv[argptr++]);
+                if (tmp < 1.0) {
+                  cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
+                  tmp = 1.3;
+                }
+//                for(int i=0;i<NDIM;i++) {
+//                  repetition[i] = atoi(argv[argptr++]);
+//                  if (repetition[i] < 1)
+//                    cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
+//                  repetition[i] = 1;
+//                }
+                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, tmp);
+              }
+              break;
+            case 'd':
+              for (int axis = 1; axis <= NDIM; axis++) {
+                int faktor = atoi(argv[argptr++]);
+                int count;
+                element ** Elements;
+                vector ** Vectors;
+                if (faktor < 1) {
+                  cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
+                  faktor = 1;
+                }
+                mol->CountAtoms((ofstream *)&cout);  // recount atoms
+                if (mol->AtomCount != 0) {  // if there is more than none
+                  count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
+                  Elements = new element *[count];
+                  Vectors = new vector *[count];
+                  j = 0;
+                  first = mol->start;
+                  while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
+                    first = first->next;
+                    Elements[j] = first->type;
+                    Vectors[j] = &first->x;
+                    j++;
+                  }
+                  if (count != j)
+                    cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+                  x.Zero();
+                  y.Zero();
+                  y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+                  for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
+                    x.AddVector(&y); // per factor one cell width further
+                    for (int k=count;k--;) { // go through every atom of the original cell
+                      first = new atom(); // create a new body
+                      first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
+                      first->x.AddVector(&x);      // translate the coordinates
+                      first->type = Elements[k];  // insert original element
+                      mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
+                    }
+                  }
+                  // free memory
+                  delete[](Elements);
+                  delete[](Vectors);
+                  // correct cell size
+                  if (axis < 0) { // if sign was negative, we have to translate everything
+                    x.Zero();
+                    x.AddVector(&y);
+                    x.Scale(-(faktor-1));
+                    mol->Translate(&x);
+                  }
+                  mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+                }
+              }
               break;
             default:   // no match? Step on
