Changeset 62d092 for src/Jobs

Timestamp:

Oct 13, 2015, 8:14:33 PM (9 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:

3690e4

Parents:

17e4fd

git-author:

Frederik Heber <heber@…> (09/08/15 22:12:20)

git-committer:

Frederik Heber <heber@…> (10/13/15 20:14:33)

Message:

Rewrote ChargeSmearer to use visitor pattern.

• this is used for going over the bspline domain.
• we need it when the support does not fit onto grid and when we need to recalculate the normalization.

Location:

src/Jobs

Files:

• ChargeSmearer.cpp (modified) (4 diffs)
• ChargeSmearer.hpp (modified) (2 diffs)

src/Jobs/ChargeSmearer.cpp

@@ -37 +37 @@
 #include "ChargeSmearer.hpp"
 
+#include "CodePatterns/Log.hpp"
 #include "CodePatterns/Singleton_impl.hpp"
 
@@ -68 +69 @@
   // reallocate memory for spline values
   delete[] vals;
-  double *vals = new double[VMG::Helper::intpow(2*nfc+1,3)];
+  vals = new double[VMG::Helper::intpow(2*nfc+1,3)];
 
   // recalculate spline values
   unsigned int c=0;
   int_val = 0.;
-  for (double dir_x=-1*nfc*meshwidth; dir_x<=nfc*meshwidth; dir_x+=meshwidth) {
-    for (double dir_y=-1*nfc*meshwidth; dir_y<=nfc*meshwidth; dir_y+=meshwidth) {
-      for (double dir_z=-1*nfc*meshwidth; dir_z<=nfc*meshwidth; dir_z+=meshwidth) {
-
+  std::pair<int, int> bounds[3];
+  for (int dim=0;dim<3;++dim) {
+    bounds[dim].first = -1*nfc;
+    bounds[dim].second = nfc;
+  }
+  for (int i=bounds[0].first; i<=bounds[0].second; ++i)
+    for (int j=bounds[1].first; j<=bounds[1].second; ++j)
+      for (int k=bounds[2].first; k<=bounds[2].second; ++k) {
+        const double dir_x = (double)i*meshwidth;
+        const double dir_y = (double)j*meshwidth;
+        const double dir_z = (double)k*meshwidth;
         // Compute distance from grid point to particle
         const double distance = std::sqrt(dir_x*dir_x+dir_y*dir_y+dir_z*dir_z);
@@ -82 +90 @@
         vals[c++] = temp_val;
         int_val += temp_val;
+        LOG(4, "DEBUG: Spline evaluated at " << distance << " is " << temp_val);
       }
-    }
-  }
 
   /// then integrate
-  int_val = 1.0 / (int_val * meshwidth * meshwidth * meshwidth);
+  int_val = 1.0 / int_val;
+
+  LOG(3, "DEBUG: Spline for smearing electronic charge distribution has norm factor of " << int_val);
+}
+
+void ChargeSmearer::visitBSplineDomain(
+    const VMG::GridIterator &_iter,
+    const visitor_t &_visitor
+    ) const
+{
+  unsigned int c = 0;
+  VMG::Index index;
+  const VMG::Index end = _iter.GetIndex() + (int)nfc;
+  for (index[0] = _iter.GetIndex()[0] - (int)nfc; index[0]<=end[0]; ++index[0])
+    for (index[1] = _iter.GetIndex()[1] - (int)nfc; index[1]<=end[1]; ++index[1])
+      for (index[2] = _iter.GetIndex()[2] - (int)nfc; index[2]<=end[2]; ++index[2]) {
+        if (index.IsInBounds(_iter.GetBegin(), _iter.GetEnd()))
+          _visitor(index, vals[c++]);
+        else
+          ++c;
+      }
+  assert( c == (unsigned int)VMG::Helper::intpow(2*nfc+1,3) );
+}
+
+static void DomainIntegrator(
+    const VMG::Index & _index,
+    const double _val,
+    double &_int_val)
+{
+  _int_val += _val;
+}
+
+static void SplineSetter(
+    const VMG::Index & _index,
+    const double _val,
+    const double _factor,
+    VMG::Grid& _grid)
+{
+  _grid(_index) += _val * _factor;
 }
 
@@ -95 +140 @@
     const double _charge) const
 {
-  /// prepare bounds
-  std::pair<int, int> bounds[3];
+  // check whether we go over whole bspline support
+  bool renormalize = false;
   for (int dim=0;dim<3;++dim) {
-    bounds[dim].first = std::max(_iter.GetBegin()[dim], _iter.GetIndex()[dim] - (int)nfc);
-    bounds[dim].second = std::min(_iter.GetEnd()[dim], _iter.GetIndex()[dim] + (int)nfc);
+    renormalize |= (_iter.GetBegin()[dim] > _iter.GetIndex()[dim] - (int)nfc);
+    renormalize |= (_iter.GetEnd()[dim] < _iter.GetIndex()[dim] + (int)nfc);
+  }
+
+  /// renormalize bspline if necessary (bounds don't fit in window)
+  double temp_int_val = 0.;
+  if (renormalize) {
+    const visitor_t integrator = boost::bind(DomainIntegrator, _1, _2, boost::ref(temp_int_val));
+    visitBSplineDomain(_iter, integrator);
+    temp_int_val = 1.0 / temp_int_val;
+  } else {
+    temp_int_val = int_val;
   }
 
   /// then transfer to grid in bounded window
-  unsigned int c = 0;
-  for (int i=bounds[0].first; i<=bounds[0].second; ++i)
-    for (int j=bounds[1].first; j<=bounds[1].second; ++j)
-      for (int k=bounds[2].first; k<=bounds[2].second; ++k)
-        _grid(i,j,k) += vals[c++] * int_val * _charge;
-  assert( c == (unsigned int)VMG::Helper::intpow(2*nfc+1,3) );
+  {
+    const visitor_t setter = boost::bind(SplineSetter, _1, _2, temp_int_val * _charge, boost::ref(_grid));
+    visitBSplineDomain(_iter, setter);
+  }
 }
 

src/Jobs/ChargeSmearer.hpp

@@ -19 +19 @@
   class Grid;
   class GridIterator;
+  class Index;
 
   namespace Particle {
@@ -109 +110 @@
 
 private:
+  //!> typedef for the visitor functor
+  typedef boost::function<void (const VMG::Index &, const double)> visitor_t;
+
+  /** Visitor pattern for going over bspline support and executing the given
+   * functor at each valid node.
+   *
+   * \param _iter iterator giving central node for support on grid
+   * \param _visitor visitor functor to execute at each node (that's still inside grid)
+   */
+  void visitBSplineDomain(
+      const VMG::GridIterator &_iter,
+      const visitor_t &_visitor
+      ) const;
+
+private:
   //!> near field radius for bspline
   unsigned int nfc;