Changeset ffe885

Timestamp:

Apr 12, 2010, 4:25:50 PM (15 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:

09898c

Parents:

4f9e47

Message:

BUGFIX: LinkedCell::GetPointsInsideSphere() did not work if center of sphere was outside of LinkedCell's domain.

• new function LinkedCell::SetClosestIndexToOutsideVector() similar to LinkedCell::SetIndexToVector(), only vector may be outside of domain. Then, closest cell is taken and minimum distance calculated.
• LinkedCell::GetPointsInsideSphere() uses LinkedCell::SetClosestIndexToOutsideVector() instead of LinkedCell::SetIndexToVector()

Location:

src

Files:

• linkedcell.cpp (modified) (2 diffs)
• linkedcell.hpp (modified) (2 diffs)

src/linkedcell.cpp

@@ -337 +337 @@
 };
 
+/** Set the index to the cell containing a given Vector *x, which is not inside the LinkedCell's domain
+ * Note that as we have to check distance from every corner of the closest cell, this function is faw more
+ * expensive and if Vector is known to be inside LinkedCell's domain, then SetIndexToVector() should be used.
+ * \param *x Vector with coordinates
+ * \return minimum squared distance of cell to given vector (if inside of domain, distance is 0)
+ */
+double LinkedCell::SetClosestIndexToOutsideVector(const Vector * const x) const
+{
+  for (int i=0;i<NDIM;i++) {
+    n[i] = (int)floor((x->x[i] - min.x[i])/RADIUS);
+    if (n[i] < 0)
+      n[i] = 0;
+    if (n[i] >= N[i])
+      n[i] = N[i]-1;
+  }
+
+  // calculate distance of cell to vector
+  double distanceSquared = 0.;
+  bool outside = true;  // flag whether x is found in- or outside of LinkedCell's domain/closest cell
+  Vector corner; // current corner of closest cell
+  Vector tester; // Vector pointing from corner to center of closest cell
+  Vector Distance;  // Vector from corner of closest cell to x
+
+  Vector center;  // center of the closest cell
+  for (int i=0;i<NDIM;i++)
+    center.x[i] = min.x[i]+((double)n[i]+.5)*RADIUS;
+
+  int c[NDIM];
+  for (c[0]=0;c[0]<=1;c[0]++)
+    for (c[1]=0; c[1]<=1;c[1]++)
+      for (c[2]=0; c[2]<=1;c[2]++) {
+        // set up corner
+        for (int i=0;i<NDIM;i++)
+          corner.x[i] = min.x[i]+RADIUS*((double)n[i]+c[i]);
+        // set up distance vector
+        Distance.CopyVector(x);
+        Distance.SubtractVector(&corner);
+        const double dist = Distance.NormSquared();
+        // check whether distance is smaller
+        if (dist< distanceSquared)
+          distanceSquared = dist;
+        // check whether distance vector goes inside or outside
+        tester.CopyVector(&center);
+        tester.SubtractVector(&corner);
+        if (tester.ScalarProduct(&Distance) < 0)
+          outside = false;
+      }
+  return (outside ? distanceSquared : 0.);
+};
 
 /** Returns a list of all TesselPoint with distance less than \a radius to \a *Center.
@@ -351 +400 @@
 
   // set index of LC to center of sphere
-  if (!SetIndexToVector(center)) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << *center << " is outside of LinkedCell's bounding box." << endl);
+  const double dist = SetClosestIndexToOutsideVector(center);
+  if (dist > radius) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << *center << " is too far away from any atom in LinkedCell's bounding box." << endl);
     return TesselList;
   }
 
   // gather all neighbours first, then look who fulfills distance criteria
-  NeighbourList = GetallNeighbours(radius);
+  NeighbourList = GetallNeighbours(radius-dist);
   if (NeighbourList != NULL) {
     for (LinkedNodes::const_iterator Runner = NeighbourList->begin(); Runner != NeighbourList->end(); Runner++) {

src/linkedcell.hpp

@@ -61 +61 @@
     bool SetIndexToNode(const TesselPoint * const Walker)const ;
     bool SetIndexToVector(const Vector * const x)const ;
+    double SetClosestIndexToOutsideVector(const Vector * const x) const;
     bool CheckBounds()const ;
     bool CheckBounds(const int relative[NDIM])const ;
@@ -67 +68 @@
     LinkedCell::LinkedNodes* GetallNeighbours(const double distance = 0) const;
     LinkedCell::LinkedNodes* GetPointsInsideSphere(const double radius, const Vector * const center) const;
-
     // not implemented yet
     bool AddNode(Vector *Walker);