Changeset 0779a9 for molecuilder

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.

Location:

molecuilder/src

Files:

• boundary.cpp (modified) (3 diffs)
• boundary.hpp (modified) (1 diff)
• builder.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;
 };
 

molecuilder/src/boundary.hpp

@@ -124 +124 @@
 double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPoints, molecule *mol);
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity);
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
 
 

molecuilder/src/builder.cpp

@@ -736 +736 @@
   int ExitFlag = 0;
   int j;
+  double volume = 0.;
   enum ConfigStatus config_present = absent;
   clock_t start,end;
@@ -987 +988 @@
               VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
               break;
+            case 'U':
+              volume = atof(argv[argptr++]);
+              cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;  
             case 'u':
               {
-                double tmp;
+                double density;
                 ExitFlag = 1;
                 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
-                tmp = atof(argv[argptr++]);
-                if (tmp < 1.0) {
+                density = atof(argv[argptr++]);
+                if (density < 1.0) {
                   cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
-                  tmp = 1.3;
+                  density = 1.3;
                 }
 //                for(int i=0;i<NDIM;i++) {
@@ -1003 +1007 @@
 //                  repetition[i] = 1;
 //                }
-                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, tmp);
+                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);
               }
               break;