Context Navigation

-              raa55d0
+              r2b7375
 };
-/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
- * We get cluster volume by Tesselation::getVolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
- * TODO: Here, we need a VolumeOfGeneralEnvelope (i.e. non-convex one)
- * \param *out output stream for debugging
- * \param *configuration needed for path to store convex envelope file
- * \param *mol molecule structure representing the cluster
- * \param *&TesselStruct Tesselation structure with triangles on return
- * \param ClusterVolume guesstimated cluster volume, if equal 0 we used Tesselation::getVolumeOfConvexEnvelope() instead.
- * \param celldensity desired average density in final cell
- */
-void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
+{
-        //Info FunctionInfo(__func__);
-  bool IsAngstroem = true;
-  double *GreatestDiameter = NULL;
-  Boundaries *BoundaryPoints = NULL;
-  class Tesselation *TesselStruct = NULL;
-  Vector BoxLengths;
-  int repetition[NDIM] = { 1, 1, 1 };
-  int TotalNoClusters = 1;
-  double totalmass = 0.;
-  double clustervolume = 0.;
-  double cellvolume = 0.;
-  // transform to PAS by Action
-  Vector MainAxis(0.,0.,1.);
-  mol->RotateToPrincipalAxisSystem(MainAxis);
-  IsAngstroem = configuration->GetIsAngstroem();
-  BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
-  GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
-  PointCloudAdaptor< molecule > cloud(mol, mol->name);
-  LinkedCell_deprecated *LCList = new LinkedCell_deprecated(cloud, 10.);
-  FindConvexBorder(mol, BoundaryPoints, TesselStruct, (const LinkedCell_deprecated *&)LCList, NULL);
-  delete (LCList);
-  delete[] BoundaryPoints;
-  // some preparations beforehand
-  if (ClusterVolume == 0)
-    clustervolume = TesselStruct->getVolumeOfConvexEnvelope(configuration->GetIsAngstroem());
-  else
-    clustervolume = ClusterVolume;
-  delete TesselStruct;
-  for (int i = 0; i < NDIM; i++)
-    TotalNoClusters *= repetition[i];
-  // sum up the atomic masses
-  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
-      totalmass += (*iter)->getType()->getMass();
+  }
-  LOG(0, "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit.");
-  LOG(0, "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
-  // solve cubic polynomial
-  LOG(1, "Solving equidistant suspension in water problem ...");
-  if (IsAngstroem)
-    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
-  else
-    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
-  LOG(1, "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
-  double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
-  LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
-  if (minimumvolume > cellvolume) {
-    ELOG(1, "the containing box already has a greater volume than the envisaged cell volume!");
-    LOG(0, "Setting Box dimensions to minimum possible, the greatest diameters.");
-    for (int i = 0; i < NDIM; i++)
-      BoxLengths[i] = GreatestDiameter[i];
-    mol->CenterEdge();
-  } else {
-    BoxLengths[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
-    BoxLengths[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[2] = minimumvolume - cellvolume;
-    double x0 = 0.;
-    double x1 = 0.;
-    double x2 = 0.;
-    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
-      LOG(0, "RESULT: The resulting spacing is: " << x0 << " .");
-    else {
-      LOG(0, "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " .");
-      x0 = x2; // sorted in ascending order
+    }
-    cellvolume = 1.;
-    for (int i = 0; i < NDIM; i++) {
-      BoxLengths[i] = repetition[i] * (x0 + GreatestDiameter[i]);
-      cellvolume *= BoxLengths[i];
+    }
-    // set new box dimensions
-    LOG(0, "Translating to box with these boundaries.");
+    {
-      RealSpaceMatrix domain;
-      for(int i =0; i<NDIM;++i)
-        domain.at(i,i) = BoxLengths[i];
-      World::getInstance().setDomain(domain);
+    }
-    mol->CenterInBox();
+  }
-  delete[] GreatestDiameter;
-  // update Box of atoms by boundary
+  {
-    RealSpaceMatrix domain;
-    for(int i =0; i<NDIM;++i)
-      domain.at(i,i) = BoxLengths[i];
-    World::getInstance().setDomain(domain);
+  }
-  LOG(0, "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
-};
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
  * \param *out output stream for debugging

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 2b7375 for src/Tesselation/boundary.cpp

Legend:

src/Tesselation/boundary.cpp

Download in other formats: