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

File:

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

src/Actions/FillAction/SuspendInMoleculeAction.cpp

@@ -129 +129 @@
   filler->CenterEdge();
 
+  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
-  std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
   double totalmass = 0.;
   const bool IsAngstroem = true;
@@ -139 +144 @@
       iter != molecules.end(); ++iter)
   {
+    // skip the filler
+    if (*iter == filler)
+      continue;
     molecule & mol = **iter;
     const double mass = calculateMass(mol);
@@ -152 +160 @@
   LOG(1, "INFO: The average density is " << setprecision(10)
       << totalmass / clustervolume << " atomicmassunit/"
-      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+      << (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;
+  }
 
   // calculate maximum solvent density
   std::vector<double> fillerdiameters(NDIM, 0.);
-  const double solventdensity =
-      calculateMass(*filler) / calculateEnvelopeVolume(*filler, fillerdiameters);
-
-  std::vector<unsigned int> counts(3, 0);
-  Vector offset(.5,.5,.5);
+  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;
 
   /// solve cubic polynomial
   double cellvolume = 0.;
   LOG(1, "Solving equidistant suspension in water problem ...");
-  cellvolume = (totalmass / solventdensity
-      - (totalmass / clustervolume)) / (params.density.get() - 1.);
+  // 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;
   LOG(1, "Cellvolume needed for a density of " << params.density.get()
       << " g/cm^3 is " << cellvolume << " angstroem^3.");
@@ -173 +188 @@
       (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) {
@@ -187 +202 @@
         + GreatestDiameter[1] * GreatestDiameter[2];
     BoxLengths[2] = minimumvolume - cellvolume;
-    double x0 = 0.;
-    double x1 = 0.;
-    double x2 = 0.;
+    std::vector<double> x(3, 0.);
     // for cubic polynomial there are either 1 or 3 unique solutions
-    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
-    }
+    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 << " .");
 
     cellvolume = 1.;
     for (size_t i = 0; i < NDIM; ++i) {
-      BoxLengths[i] = x0 + GreatestDiameter[i];
+      BoxLengths[i] = x[i] + GreatestDiameter[i];
       cellvolume *= BoxLengths[i];
     }
-
-    // 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();
+  }
+
+  // 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);
   }
 
@@ -231 +244 @@
       params.RandMoleculeDisplacement.get(),
       params.DoRotate.get());
+  Vector offset(.5,.5,.5);
   filler_preparator.addCubeMesh(
       counts,