Changeset 1d2d01e for src/LinkedCell

Timestamp:

Apr 17, 2013, 6:56:06 PM (12 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:

79d0b9

Parents:

0773bd

git-author:

Frederik Heber <heber@…> (03/18/13 18:39:07)

git-committer:

Frederik Heber <heber@…> (04/17/13 18:56:06)

Message:

Fixing minimum sensible length for LinkedCell_Controller.

• we did use the minimum sensible length for getBestModel() but not for the actual creation.
• extracted code into function getSensibleMinimumLength().
• fixed (and documented) formula for calculating min_distance.
• we also can now also deal with systems that are mostly 2d or 1d.

Location:

src/LinkedCell

Files:

• LinkedCell_Controller.cpp (modified) (5 diffs)
• LinkedCell_Controller.hpp (modified) (1 diff)

src/LinkedCell/LinkedCell_Controller.cpp

@@ -116 +116 @@
 }
 
-/** Internal function to decide whether a suitable model is present or not.
- *
- * Here, the heuristic for deciding whether a new linked cell structure has to
- * be constructed or not is implemented. The current heuristic is as follows:
- * -# the best model should have at least half the desired length (such
- *    that most we have to look two neighbor shells wide and not one).
- * -# the best model should have at most twice the desired length but
- *    no less than 1 angstroem.
- *
- * \note Dealing out a pointer is here (hopefully) safe because the function is
- * internal and we - inside this class - know what we are doing.
- *
- * @param distance edge length of the requested linked cell structure
- * @return NULL - there is no fitting LinkedCell_Model, else - pointer to instance
- */
-const LinkedCell_Model *LinkedCell_Controller::getBestModel(double distance) const
+static size_t checkLengthsLess(
+    const Vector &Lengths,
+    const double min_distance)
+{
+  size_t counter=0;
+  for (size_t i=0;i<NDIM;++i) {
+    if (Lengths[i] < min_distance)
+      ++counter;
+  }
+  return counter;
+}
+
+/** Given some given \a distance return it such that it is larger equal than some minimum
+ * sensible value with respect to the current domain.
+ *
+ * The sensible minimum distance \f$\Delta\f$ is computed as follows:
+ * \f$ \frac {l_1} \Delta \frac {l_2} \Delta \frac {l_3} \Delta \leq M \f$, where M is the
+ * maximum number of linked cells. Where some additional thought must be given
+ * to rounding issues. Hence, we calculate \f$ (\frac {l_i} d)^{\frac 1 d}+.5\f$,
+ * where d is the assumed dimensionality of the domain: d=3 for cubic, d=2 for
+ * a planar, d=1 for a bar-like system.
+ *
+ * \param distance distance compared with minimum value
+ * \return minimum sensible value or \a distance if larger
+ */
+double LinkedCell_Controller::getMinimumSensibleLength(double distance) const
 {
   // first check the box with respect to the given distance
@@ -141 +151 @@
   UnitMatrix *= M;
   Vector Lengths = UnitMatrix.transformToEigenbasis();
-  double min_distance = std::numeric_limits<double>::max();
   double max_length = 0.;
-  for (size_t i=0;i<NDIM;++i) {
-    min_distance = std::min(Lengths[i], min_distance);
-    max_length = std::max(Lengths[i], max_length);
-  }
-  min_distance *= 1e-3; // at most 1000 cells per axis
+  double volume = 1.;
+  double min_distance = lower_threshold;
+  for (size_t dim=NDIM; dim>0; --dim) {
+    const double fraction = pow(MAX_LINKEDCELLNODES+1., 1./(double)dim);
+    for (size_t i=0;i<NDIM;++i) {
+      max_length = std::max(Lengths[i], max_length);
+      // the length dropping out for dimensionality less than 3 must appear
+      // just as factor of 1 not less.
+      if (Lengths[i] < min_distance)
+        volume *= pow(min_distance, 1./(double)dim)+.5;
+      else
+        volume *= pow(Lengths[i], 1./(double)dim)+.5; //.5 for rounding issues
+    }
+    min_distance = volume/fraction;
+    // assuming a (cubic-like) 3d system
+    if (checkLengthsLess(Lengths, min_distance) == 0)
+      break;
+  }
   if (distance < min_distance) {
     LOG(1, "INFO: Setting LC distance from given " << distance
@@ -153 +175 @@
     distance = min_distance;
   }
-
+  return distance;
+}
+
+/** Internal function to decide whether a suitable model is present or not.
+ *
+ * Here, the heuristic for deciding whether a new linked cell structure has to
+ * be constructed or not is implemented. The current heuristic is as follows:
+ * -# the best model should have at least half the desired length (such
+ *    that most we have to look two neighbor shells wide and not one).
+ * -# the best model should have at most twice the desired length but
+ *    no less than 1 angstroem.
+ *
+ * \note Dealing out a pointer is here (hopefully) safe because the function is
+ * internal and we - inside this class - know what we are doing.
+ *
+ * @param distance edge length of the requested linked cell structure
+ * @return NULL - there is no fitting LinkedCell_Model, else - pointer to instance
+ */
+const LinkedCell_Model *LinkedCell_Controller::getBestModel(double distance) const
+{
   /// Bound distance to be within [lower_threshold, upper_threshold).
   /// Note that we need to stay away from upper boundary a bit,
@@ -159 +200 @@
   if (distance < lower_threshold)
     distance = lower_threshold;
-  if (distance > max_length)
-    distance = max_length - std::numeric_limits<double>::round_error();
+  if (distance > upper_threshold)
+    distance = upper_threshold - std::numeric_limits<double>::round_error();
 
   /// Look for all models within [0.5 distance, 2. distance).
@@ -210 +251 @@
 LinkedCell_View LinkedCell_Controller::getView(const double distance, IPointCloud &set)
 {
+  // distance should be a given minimum length dependent on domain at least
+  const double ConstraintDistance = getMinimumSensibleLength(distance);
+
   /// Look for best instance.
-  const LinkedCell_Model * const LCModel_best = getBestModel(distance);
+  const LinkedCell_Model * const LCModel_best = getBestModel(ConstraintDistance);
 
   /// Construct new instance if none found,
   if (LCModel_best == NULL) {
-    LinkedCell_Model * const LCModel_new = new LinkedCell_Model(distance, domain);
+    LinkedCell_Model * const LCModel_new = new LinkedCell_Model(ConstraintDistance, domain);
     LCModel_new->insertPointCloud(set);
-    insertNewModel(distance, LCModel_new);
+    insertNewModel(ConstraintDistance, LCModel_new);
     LinkedCell_View view(*LCModel_new);
     return view;

src/LinkedCell/LinkedCell_Controller.hpp

@@ -61 +61 @@
     void updateModels();
 
+    double getMinimumSensibleLength(double distance) const;
+
   private:
     //!> typedef for the (edge length, model) map