Changeset 2440ce for src/Actions/FillAction/SuspendInMoleculeAction.cpp

Timestamp:

Jan 14, 2015, 9:03:24 PM (10 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

cad383

Parents:

601ef8

git-author:

Frederik Heber <heber@…> (01/14/15 20:31:24)

git-committer:

Frederik Heber <heber@…> (01/14/15 21:03:24)

Message:

FIX: Changed SuspendInMoleculeAction to catch segfault when rho=1 was given.

• however, the action is still not tested to work.
• TESTFIX: suspend-in-water gets parameter density via "density" not via "suspend-in-water", set desired density to faulting 1 for the moment.
• TESTFIX: added water molecule to test.conf, to have at least two molecules.

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,