Changes in src/Actions/FillAction/SuspendInMoleculeAction.cpp [2440ce:aa55d0]

File:

• src/Actions/FillAction/SuspendInMoleculeAction.cpp (modified) (6 diffs)

src/Actions/FillAction/SuspendInMoleculeAction.cpp

@@ -129 +129 @@
   filler->CenterEdge();
 
+  /// first we need to calculate some volumes and masses
   std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
-  if (molecules.size() < 2) {
-    STATUS("There must be at least two molecules: filler and to be suspended.");
-    return Action::failure;
-  }
-
-  /// first we need to calculate some volumes and masses
   double totalmass = 0.;
   const bool IsAngstroem = true;
@@ -144 +139 @@
       iter != molecules.end(); ++iter)
   {
-    // skip the filler
-    if (*iter == filler)
-      continue;
     molecule & mol = **iter;
     const double mass = calculateMass(mol);
@@ -160 +152 @@
   LOG(1, "INFO: The average density is " << setprecision(10)
       << totalmass / clustervolume << " atomicmassunit/"
-      << (IsAngstroem ? " angstrom" : " atomiclength") << "^3.");
-  if ( ((totalmass / clustervolume < 1.) && (params.density.get() > 1.))
-      || ((totalmass / clustervolume > 1.) && (params.density.get() < 1.))) {
-    STATUS("Desired and present molecular densities must both be either in [0,1) or in (1, inf).");
-    return Action::failure;
-  }
+      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
 
   // calculate maximum solvent density
   std::vector<double> fillerdiameters(NDIM, 0.);
-  const double fillervolume = calculateEnvelopeVolume(*filler, fillerdiameters);
-  const double fillermass = calculateMass(*filler);
-  LOG(1, "INFO: The filler's mass is " << setprecision(10)
-      << fillermass << " atomicmassunit, and it's volume is "
-      << fillervolume << (IsAngstroem ? " angstrom" : " atomiclength") << "^3.");
-  const double solventdensity = fillermass / fillervolume;
+  const double solventdensity =
+      calculateMass(*filler) / calculateEnvelopeVolume(*filler, fillerdiameters);
+
+  std::vector<unsigned int> counts(3, 0);
+  Vector offset(.5,.5,.5);
 
   /// solve cubic polynomial
   double cellvolume = 0.;
   LOG(1, "Solving equidistant suspension in water problem ...");
-  // s = solvent, f = filler, 0 = initial molecules/cluster
-  // v_s = v_0 + v_f, m_s = m_0 + rho_f * v_f --> rho_s = m_s/v_s ==> v_f = (m_0 - rho_s * v_o) / (rho_s - rho_f)
-  cellvolume = (totalmass - params.density.get() * clustervolume) / (params.density.get() - 1.) + clustervolume;
+  cellvolume = (totalmass / solventdensity
+      - (totalmass / clustervolume)) / (params.density.get() - 1.);
   LOG(1, "Cellvolume needed for a density of " << params.density.get()
       << " g/cm^3 is " << cellvolume << " angstroem^3.");
@@ -188 +173 @@
       (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
   LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is "
-      << minimumvolume << " angstrom^3.");
+      << minimumvolume << "angstrom^3.");
 
   if (minimumvolume > cellvolume) {
@@ -202 +187 @@
         + GreatestDiameter[1] * GreatestDiameter[2];
     BoxLengths[2] = minimumvolume - cellvolume;
-    std::vector<double> x(3, 0.);
+    double x0 = 0.;
+    double x1 = 0.;
+    double x2 = 0.;
     // for cubic polynomial there are either 1 or 3 unique solutions
-    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x[0], &x[1], &x[2]) == 1) {
-      x[1] = x[0];
-      x[2] = x[0];
-    } else {
-      std::swap(x[0], x[2]); // sorted in ascending order
-    }
-    LOG(0, "RESULT: The resulting spacing is: " << x << " .");
+    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1)
+      LOG(0, "RESULT: The resulting spacing is: " << x0 << " .");
+    else {
+      LOG(0, "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " .");
+      std::swap(x0, x2); // sorted in ascending order
+    }
 
     cellvolume = 1.;
     for (size_t i = 0; i < NDIM; ++i) {
-      BoxLengths[i] = x[i] + GreatestDiameter[i];
+      BoxLengths[i] = x0 + GreatestDiameter[i];
       cellvolume *= BoxLengths[i];
     }
-  }
-
-  // TODO: Determine counts from resulting mass correctly (hard problem due to integers)
-  std::vector<unsigned int> counts(3, 0);
-  const unsigned int totalcounts = round(params.density.get() * cellvolume - totalmass) / fillermass;
-  if (totalcounts > 0) {
-    counts[0] = ceil(BoxLengths[0]/3.1);
-    counts[1] = ceil(BoxLengths[1]/3.1);
-    counts[2] = ceil(BoxLengths[2]/3.1);
+
+    // set new box dimensions
+    LOG(0, "Translating to box with these boundaries.");
+    {
+      RealSpaceMatrix domain;
+      for(size_t i =0; i<NDIM;++i)
+        domain.at(i,i) = BoxLengths[i];
+      World::getInstance().setDomain(domain);
+    }
+//    mol->CenterInBox();
   }
 
@@ -244 +231 @@
       params.RandMoleculeDisplacement.get(),
       params.DoRotate.get());
-  Vector offset(.5,.5,.5);
   filler_preparator.addCubeMesh(
       counts,