Changeset 160ad7 for src/Fragmentation

Timestamp:

Sep 14, 2016, 6:42:53 PM (8 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, 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_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable

Children:

ec1aae

Parents:

b839fb

git-author:

Frederik Heber <heber@…> (03/10/16 15:34:11)

git-committer:

Frederik Heber <heber@…> (09/14/16 18:42:53)

Message:

ExportGraph_ToJobs::operator() now calculates a specific level for each fragment.

• this is based on a set maximum mesh width that would still be acceptable.
• added unit test on introduced static helper functions for ExportGraph_ToJobs.
• added assert to make sure desired_level is never negative.
• we check that conditions are met for resulting level inside getGridExtentsForGivenBoundingBox().
• added extensive explanation.

Location:

src/Fragmentation/Exporters

Files:

• ExportGraph_ToJobs.cpp (modified) (6 diffs)
• ExportGraph_ToJobs.hpp (modified) (4 diffs)
• unittests/ExportGraph_ToJobsUnitTest.cpp (added)
• unittests/ExportGraph_ToJobsUnitTest.hpp (added)
• unittests/Makefile.am (modified) (4 diffs)

src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

@@ -42 +42 @@
 
 #include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include "CodePatterns/Log.hpp"
 
 #include "Box.hpp"
@@ -56 +60 @@
 #include "World.hpp"
 
+const double ExportGraph_ToJobs::log_two(log(2.));
+
 ExportGraph_ToJobs::ExportGraph_ToJobs(
     const Graph &_graph,
@@ -62 +68 @@
     const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
     ExportGraph(_graph, _treatment, _saturation,_globalsaturationpositions),
-    level(5)
+    level(5),
+    max_meshwidth(0.),
+    minimum_empty_boundary(5./AtomicLengthToAngstroem)
 {}
 
@@ -80 +88 @@
 }
 
+SamplingGridProperties ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox(
+    const std::pair<Vector, Vector> &_minmax,
+    const SamplingGridProperties &_domain,
+    const double & _minimum_empty_boundary
+    )
+{
+  // return value
+  SamplingGridProperties fragment_window;
+
+  /// determine center of fragment
+  const Vector center = .5*(_minmax.first + _minmax.second);
+  LOG(4, "DEBUG: center of fragment is at " << center);
+
+  /// associate center to its containing grid cell (defined by boundary points)
+  Vector lower_center;
+  Vector higher_center;
+  bool equal_components[NDIM] = { false, false, false };
+  for (size_t i=0;i<NDIM;++i) {
+    lower_center[i] = _domain.getNearestLowerGridPoint(center[i], i);
+    higher_center[i] = _domain.getNearestHigherGridPoint(center[i], i);
+    equal_components[i] =
+        fabs(lower_center[i] - higher_center[i]) < std::numeric_limits<double>::epsilon()*1e4;
+  }
+  LOG(5, "DEBUG: lower_center is " << lower_center << ", higher_center is " << higher_center);
+
+  // fashion min max into cubic box extents of at least extent plus empty
+  // boundary and at most three times the extent
+  const Vector extent = _minmax.second - _minmax.first;
+  double greatest_extent = extent[extent.GreatestComponent()];
+  if (greatest_extent > _minimum_empty_boundary)
+    greatest_extent *= 3.;
+  else
+    greatest_extent += 2.* _minimum_empty_boundary;
+  const Vector total(
+      _domain.getTotalLengthPerAxis(0),
+      _domain.getTotalLengthPerAxis(1),
+      _domain.getTotalLengthPerAxis(2)
+  );
+  const double greatest_total = total[total.GreatestComponent()];
+  greatest_extent = std::min(greatest_extent, greatest_total);
+  LOG(5, "DEBUG: extent of fragment is " << extent << ", greatest_extent is " << greatest_extent);
+
+  /// increase levels around this center to find the matching window
+  const double delta[NDIM] = {
+      _domain.getDeltaPerAxis(0),
+      _domain.getDeltaPerAxis(1),
+      _domain.getDeltaPerAxis(2)
+  };
+  LOG(6, "DEBUG: delta is " << Vector(delta));
+  const double greatest_delta = std::max(delta[0], std::max(delta[1], delta[2]));
+
+  fragment_window.level = (int)ceil(log(greatest_extent/greatest_delta+1)/log_two);
+  const size_t half_fragment_gridpoints = 1 << (fragment_window.level-1);
+  const size_t domain_gridpoints = _domain.getGridPointsPerAxis();
+  int begin_index[NDIM];
+  int end_index[NDIM];
+  int begin_steps[NDIM] = { 0, 0, 0 };
+  int end_steps[NDIM] = { 0, 0, 0 };
+  for (size_t i=0;i<NDIM;++i) {
+    begin_index[i] = round(lower_center[i]/delta[i]) - (half_fragment_gridpoints-1);
+    end_index[i] = round(higher_center[i]/delta[i]) + (half_fragment_gridpoints+(unsigned int)(equal_components[i]));
+    if (begin_index[i] < 0)
+      begin_steps[i] = -begin_index[i];
+    if (end_index[i] > (int)domain_gridpoints)
+      end_steps[i] = end_index[i] - domain_gridpoints;
+    if ((begin_steps[i]+end_steps[i]+end_index[i]-begin_index[i] > (int)domain_gridpoints)
+        || ((begin_steps[i] != 0) && (end_steps[i] != 0))) {
+      begin_index[i] = 0;
+      end_index[i] = domain_gridpoints;
+      begin_steps[i] = 0;
+      end_steps[i] = 0;
+    }
+  }
+  for (size_t i=0;i<NDIM;++i) {
+    fragment_window.begin[i] = (begin_index[i]+begin_steps[i]-end_steps[i]) * delta[i];
+    fragment_window.end[i] = (end_index[i]-end_steps[i]+begin_steps[i]) * delta[i];
+  }
+    LOG(6, "DEBUG: fragment begin is " << Vector(fragment_window.begin)
+        << ", fragment end is " << Vector(fragment_window.end));
+
+    /// check whether window is large enough, if not yet continue
+#ifndef NDEBUG
+    for (size_t i=0;i<NDIM;++i) {
+      const double window_length = fragment_window.end[i] - fragment_window.begin[i];
+      ASSERT (((greatest_total != greatest_extent) && (greatest_extent+delta[i] <= window_length))
+          || ((greatest_total == greatest_extent) && (greatest_extent <= window_length)),
+          "ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox() - level "
+          +toString(fragment_window.level)+" is insufficient to place fragment inside box "
+          +toString(_minmax.first)+" <-> "+toString(_minmax.second));
+    }
+#endif
+#ifndef NDEBUG
+  for (size_t i=0;i<NDIM;++i) {
+    ASSERT( fragment_window.begin[i] >= _domain.begin[i],
+        "ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox() - fragment begin "
+        +toString(fragment_window.begin[i])+" is smaller than that of domain "
+        +toString(_domain.begin[i]));
+    ASSERT( fragment_window.end[i] <= _domain.end[i],
+        "ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox() - fragment end "
+        +toString(fragment_window.end[i])+" is smaller than that of domain "
+        +toString(_domain.end[i]));
+  }
+#endif
+  if (DoLog(6))
+    for (size_t i=0;i<NDIM;++i)
+      LOG(6, "DEBUG: We have " << (fragment_window.end[i] - fragment_window.begin[i])/delta[i]
+          << " gridpoints on axis " << i);
+
+  return fragment_window;
+}
+
+void ExportGraph_ToJobs::setAcceptableFragmentLevel(
+    SamplingGridProperties &_grid,
+    const double &_max_meshwidth
+    )
+{
+  double max_domain_meshwidth = _max_meshwidth;
+  for (size_t i=0;i<NDIM;++i)
+    max_domain_meshwidth = std::max(max_domain_meshwidth, _grid.getDeltaPerAxis(i));
+  // find power of 2 that's just greater than that ratio of given over desired mesh width
+  const int desired_level = ceil(log(max_domain_meshwidth / _max_meshwidth)/log_two);
+  // we may never make level smaller (needs to be as large as domain's)
+  ASSERT( desired_level >= 0,
+      "ExportGraph_ToJobs::setAcceptableFragmentLevel() - should never get negative extra level.");
+  if (desired_level > 0)
+    _grid.level +=  desired_level;
+  LOG(4, "DEBUG: Fragment requires " << desired_level
+      << " additional grid levels to reach at least " << _max_meshwidth
+      << " from " << max_domain_meshwidth);
+}
+
 bool ExportGraph_ToJobs::operator()()
 {
@@ -99 +238 @@
     const KeySet &set = CurrentFragment->getKeySet();
     LOG(2, "INFO: Creating bond fragment job for set " << set << ".");
+
+    SamplingGridProperties fragment_window(grid);
+    if (max_meshwidth != 0.) {
+      // gather extent of the fragment in atomic length(!)
+      std::pair<Vector, Vector> minmax = CurrentFragment->getRoughBoundingBox();
+      minmax.first *= 1./AtomicLengthToAngstroem;
+      minmax.second *= 1./AtomicLengthToAngstroem;
+
+      /** we need to find the begin and end points of the smaller grid
+       * otherwise we would calculate the whole domain with an incrased grid level with just
+       * a small window (which is only used for sampling and storing; vmg uses full grid).
+       * Hence, we need to make the grid smaller and such that it is still in a power of two
+       * and coinciding with the grid points of the global grid.
+       */
+      fragment_window = getGridExtentsForGivenBoundingBox(minmax, grid, minimum_empty_boundary);
+      LOG(4, "DEBUG: Fragment " << CurrentFragment->getKeySet() << " has window from  "
+          << Vector(fragment_window.begin) << " to " << Vector(fragment_window.end)
+          << " with total level portion of " << fragment_window.level );
+
+      // next we need to check how much we need to increase from the grid level for
+      // the total domain to achieve at least maximum mesh width
+      setAcceptableFragmentLevel(fragment_window, max_meshwidth/AtomicLengthToAngstroem);
+    }
+
     // store config in stream
     {
@@ -107 +270 @@
       FragmentJob::ptr testJob(
 #ifdef HAVE_JOBMARKET
-          new MPQCJob(JobId::IllegalJob, output.str(), grid.begin, grid.end, level)
+          new MPQCJob(
+              JobId::IllegalJob,
+              output.str(),
+              fragment_window.begin, fragment_window.end, fragment_window.level)
 #else
           new MPQCCommandJob(output.str(), JobId::IllegalJob)

src/Fragmentation/Exporters/ExportGraph_ToJobs.hpp

@@ -20 +20 @@
 #include "Fragmentation/Exporters/ExportGraph.hpp"
 #include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"
+#include "LinearAlgebra/Vector.hpp"
+
+class ExportGraph_ToJobsTest;
 
 /** ExportGraph_ToJobs implements an ExportGraph by sending the created
@@ -27 +30 @@
 class ExportGraph_ToJobs : public ExportGraph
 {
+  //!> grant unit test access
+  friend class ExportGraph_ToJobsTest;
 public:
   /** Constructor for ExportGraph_ToJobs.
@@ -47 +52 @@
         bool operator()();
 
-  /** Sets the level for the sampling of the density.
+        /** Sets the level for the sampling of the density.
+         *
+         * \param _level level to set
+         */
+  void setLevel(const size_t _level) { level = _level; }
+
+  /** Sets how far apart discrete points may be at most per axis.
    *
-   * \param _level level to set
+   * \param _max_meshwidth maximum allowed mesh width.
    */
-  void setLevel(const size_t _level) { level = _level; }
+  void setMaximumMeshWidth(const double _max_meshwidth) { max_meshwidth = _max_meshwidth; }
 
   /** Helper to get the global domain grid from the current simulation box.
@@ -63 +74 @@
 
 private:
+
+  /** Helper function to get the discrete extent of the grid that
+   * captures the whole of the fragment inside \a _minmax, and some
+   * additional boundary \a _minimum_empty_boundary for a given
+   * domain \a _domain with some grid level. If the fragment is larger than
+   * the empty boundary, we make the small fragment grid at most three times
+   * the extent of the fragment.
+   *
+   * We need to maintain the following properties:
+   * -# the fragment grid's begin and end need to reside (exactly) on grid points
+   *    of the global \a _domain
+   * -# the length of the fragment grid in \a _domain's deltas needs to be a
+   *    power of 2.
+   *
+   * In order to achieve these, we use the center of the fragment obtained from
+   * its extensions in \a _minmax and convert it to the next lower and upper grid
+   * points on \a _domain. For increasing powers of 2 we check the extent of the
+   * gridpoints by going half of the gridpoints back and forth relative to these
+   * center grid points. If their length relative to \a _domain's delta is
+   * sufficient to capture the desired extent of fragment, namely \a _minmax
+   * including \a _minimum_empty_boundary, we are done. We need to do this
+   * iteratively, as the fragment grid may exceed \a _domain's begin or end
+   * and shifting is required.
+   *
+   * \param _minmax minimum and maximum components of fragment's bounding box
+   * \param _domain grid with begin and end components and grid level
+   * \param _minimum_empty_boundary additional empty boundary around fragment
+   * \return grid with begin and end points and and grid level to have same spacing as domain
+   */
+  static SamplingGridProperties getGridExtentsForGivenBoundingBox(
+      const std::pair<Vector, Vector> &_minmax,
+      const SamplingGridProperties &_domain,
+      const double & _minimum_empty_boundary);
+
+  /** Helper function to calculate the required grid level for a fragment \a _grid
+   * such that is discretization mesh width is below the acceptable \a _max_meshwidth.
+   *
+   * \param _grid grid properties of the fragment, sets its level on return
+   * \param _max_meshwidth uppermost acceptable mesh width
+   */
+  static void setAcceptableFragmentLevel(
+      SamplingGridProperties &_grid,
+      const double &_max_meshwidth);
+
+  //>! log of two
+  static const double log_two;
+
+private:
   //!> resolution of sampled electron density as \f$2^{\mathrm{level}}\f$
   size_t level;
+  //!> maximum allowed mesh width, i.e. discrete points may be at most that far part per axis
+  double max_meshwidth;
+  //!> desired minimum empty boundary around the fragment
+  const double minimum_empty_boundary;
 };
 

src/Fragmentation/Exporters/unittests/Makefile.am

@@ -3 +3 @@
 
 FRAGMENTATIONEXPORTERSSOURCES = \
+        ../Fragmentation/Exporters/unittests/ExportGraph_ToJobsUnitTest.cpp \
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.cpp \
         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp \
@@ -9 +10 @@
 
 FRAGMENTATIONEXPORTERSTESTSHEADERS = \
+        ../Fragmentation/Exporters/unittests/ExportGraph_ToJobsUnitTest.hpp \
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.hpp \
         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp \
@@ -16 +18 @@
 
 FRAGMENTATIONEXPORTERSTESTS = \
+        ExportGraph_ToJobsUnitTest \
         HydrogenPoolUnitTest \
         SaturatedFragmentUnitTest \
@@ -31 +34 @@
         ${CodePatterns_LIBS} \
         $(BOOST_LIB)
+
+ExportGraph_ToJobsUnitTest_SOURCES = \
+        ../Fragmentation/Exporters/unittests/ExportGraph_ToJobsUnitTest.cpp \
+        ../Fragmentation/Exporters/unittests/ExportGraph_ToJobsUnitTest.hpp
+ExportGraph_ToJobsUnitTest_LDADD = \
+        ../libMolecuilderUI.la \
+        ../libMolecuilderFragmentationSetValues.la \
+        ${ALLLIBS}
 
 HydrogenPoolUnitTest_SOURCES = \