Changeset 0779a9 for molecuilder/src/boundary.cpp

Timestamp:

Jun 13, 2008, 2:26:50 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

1f9f1b

Parents:

3f0c46

Message:

some mor changes to fine-tune the cluster generation.

Sadly, VolumeOfConvexEnvelope() is not robust/stable, with Ratio1.4-1/1800K/compundH.*.xyz no usable convex envelope was found but points were connected right through the cluster.

File:

• molecuilder/src/boundary.cpp (modified) (3 diffs)

molecuilder/src/boundary.cpp

@@ -561 +561 @@
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
+ * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
  * \param celldensity desired average density in final cell
  */
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity)
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity)
 {
   // transform to PAS
@@ -571 +572 @@
   bool IsAngstroem = configuration->GetIsAngstroem();
   Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
-  double clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
+  double clustervolume;
+  if (ClusterVolume == 0)
+    clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
+  else 
+    clustervolume = ClusterVolume;
   double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
   vector BoxLengths;
@@ -601 +606 @@
   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;
-  
-  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]);
-  BoxLengths.x[2] = minimumvolume - cellvolume;
-  double x0 = 0.,x1 = 0.,x2 = 0.;
-  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 {
-    *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
-    x0 = x2;  // sorted in ascending order
-  }
-
-  cellvolume = 1;
-  for(int i=0;i<NDIM;i++) {
-    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);
+    cout << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
+    for(int i=0;i<NDIM;i++)
+      BoxLengths.x[i] = GreatestDiameter[i];
+    mol->CenterEdge(out, &BoxLengths);
+  } else {
+    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]);
+    BoxLengths.x[2] = minimumvolume - cellvolume;
+    double x0 = 0.,x1 = 0.,x2 = 0.;
+    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 {
+      *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
+      x0 = x2;  // sorted in ascending order
+    }
+  
+    cellvolume = 1;
+    for(int i=0;i<NDIM;i++) {
+      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+      cellvolume *= BoxLengths.x[i];
+    }
+  
+    // 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);
-  
+  *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;
 };