Changeset 826e8c for src/Graph

Timestamp:

Feb 22, 2012, 11:29:04 AM (13 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:

153985

Parents:

cbcda6

git-author:

Frederik Heber <heber@…> (02/06/12 17:46:56)

git-committer:

Frederik Heber <heber@…> (02/22/12 11:29:04)

Message:

Introduced new LinkedCell into BondGraph.

• new convenience functions BondGraph::getDomain(), ::getLinkedCell(), ::getTime() to have World dependency not in header file (circular include).
• new function BondGraph::getMaxPossibleBondDistance() where one element is given.
• dropped BondGraph::CreateAdjacency(LC_deprecated) as functionality is better placed directly into templated version directly.
• rewrote templated BondGraph::CreateAdjacency that now makes use of LinkedCell_View::getAllNeighbours().
• FIX: BondGraph::CreateAdjacency() created bonds for atom's fathers, not for atoms.
• TESTFIX: New LinkedCell construct in BondGraph::CreateAdjacency changes order of ids in adjacency file:
  • ids are now in order (unlike to before), hence we simply use the new file in the regression test.
  • this is actually a fault of the test, not of the code as order of ids is not guaranteed in the file anyway.
• TESTFIX: Removed XFAILs from regression test Graph/SubgraphDissection- BoundaryCondition.
• TESTFIX: CountBondsUnitTest use two molecules that are on top of each other: This causes mess with changed LinkedCell BondGraph recognition, and we even actually need just one molecule for the testing.

Location:

src/Graph

Files:

• BondGraph.cpp (modified) (2 diffs)
• BondGraph.hpp (modified) (7 diffs)

src/Graph/BondGraph.cpp

@@ -98 +98 @@
     int secondZ) const
 {
+  double return_length;
+  if ((firstZ < 0) || (firstZ >= (int)BondLengthMatrix->Matrix[0].size()))
+    return -1.;
+  if ((secondZ < 0) || (secondZ >= (int)BondLengthMatrix->Matrix[0][firstZ].size()))
+    return -1.;
   if (BondLengthMatrix == NULL) {
-    std::cout << "-1." << std::endl;
-    return( -1. );
+    return_length = -1.;
   } else {
-    std::cout << BondLengthMatrix->Matrix[0][firstZ][secondZ] << std::endl;
-    return (BondLengthMatrix->Matrix[0][firstZ][secondZ]);
-  }
+    return_length = BondLengthMatrix->Matrix[0][firstZ][secondZ];
+  }
+  LOG(4, "INFO: Request for length between " << firstZ << " and "
+      << secondZ << ": " << return_length << ".");
+  return return_length;
 }
 
@@ -168 +174 @@
 }
 
-void BondGraph::CreateAdjacency(LinkedCell_deprecated &LC) const
-{
-  atom *Walker = NULL;
-  atom *OtherWalker = NULL;
-  int n[NDIM];
-  Box &domain = World::getInstance().getDomain();
-  size_t CurrentTime = WorldTime::getTime();
-
-  unsigned int BondCount = 0;
-  // 3a. go through every cell
-  LOG(3, "INFO: Celling ... ");
-  for (LC.n[0] = 0; LC.n[0] < LC.N[0]; LC.n[0]++)
-    for (LC.n[1] = 0; LC.n[1] < LC.N[1]; LC.n[1]++)
-      for (LC.n[2] = 0; LC.n[2] < LC.N[2]; LC.n[2]++) {
-        const TesselPointSTLList *List = LC.GetCurrentCell();
-        LOG(2, "Current cell is " << LC.n[0] << ", " << LC.n[1] << ", " << LC.n[2] << " with No. " << LC.index << " containing " << List->size() << " points.");
-        if (List != NULL) {
-          for (TesselPointSTLList::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            Walker = dynamic_cast<atom*>(*Runner);
-            ASSERT(Walker != NULL,
-                "BondGraph::CreateAdjacency() - Tesselpoint that was not an atom retrieved from LinkedNode");
-            LOG(2, "INFO: Current Atom is " << *Walker << ".");
-            // 3c. check for possible bond between each atom in this and every one in the 27 cells
-            for (n[0] = -1; n[0] <= 1; n[0]++)
-              for (n[1] = -1; n[1] <= 1; n[1]++)
-                for (n[2] = -1; n[2] <= 1; n[2]++) {
-                  const TesselPointSTLList *OtherList = LC.GetRelativeToCurrentCell(n);
-                  if (OtherList != NULL) {
-                    LOG(3, "INFO: Current relative cell is " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << LC.index << " containing " << OtherList->size() << " points.");
-                    for (TesselPointSTLList::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
-                      if ((*OtherRunner) > Walker) {  // just to not add bonds from both sides
-                        OtherWalker = dynamic_cast<atom*>(*OtherRunner);
-                        ASSERT(OtherWalker != NULL,
-                            "BondGraph::CreateAdjacency() - TesselPoint that was not an atom retrieved from LinkedNode");
-                        LOG(3, "INFO: Current other atom is " << *OtherWalker << ".");
-                        if (OtherWalker->father > Walker->father ) { // just to not add bonds from both sides
-                          const range<double> MinMaxDistanceSquared(
-                              getMinMaxDistanceSquared(Walker, OtherWalker));
-                          const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
-                          LOG(3, "INFO: Checking squared distance " << distance << " against typical bond length of " << MinMaxDistanceSquared << ".");
-                          const bool status = MinMaxDistanceSquared.isInRange(distance);
-                          if (status) { // create bond if distance is smaller
-                            LOG(1, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
-                            //const bond * Binder =
-                              Walker->father->addBond(CurrentTime, OtherWalker->father);
-                            BondCount++;
-                          } else {
-                            LOG(2, "REJECT: Squared distance "
-                                << distance << " is out of squared covalent bounds "
-                                << MinMaxDistanceSquared << ".");
-                          }
-                        } else {
-                          LOG(5, "REJECT: Not Adding: Wrong order of father's.");
-                        }
-                      } else {
-                        LOG(5, "REJECT: Not Adding: Wrong order.");
-                      }
-                    }
-                  }
-                }
-          }
-        }
-      }
-    LOG(1, "I detected " << BondCount << " bonds in the molecule.");
+LinkedCell::LinkedCell_View BondGraph::getLinkedCell(const double max_distance) const
+{
+  return World::getInstance().getLinkedCell(max_distance);
+}
+
+Box &BondGraph::getDomain() const
+{
+  return World::getInstance().getDomain();
+}
+
+unsigned int BondGraph::getTime() const
+{
+  return WorldTime::getTime();
 }
 

src/Graph/BondGraph.hpp

@@ -24 +24 @@
 #include "Atom/AtomSet.hpp"
 #include "Bond/bond.hpp"
+#include "Box.hpp"
 #include "CodePatterns/Assert.hpp"
 #include "CodePatterns/Log.hpp"
@@ -29 +30 @@
 #include "Element/element.hpp"
 #include "Fragmentation/MatrixContainer.hpp"
-#include "LinkedCell/linkedcell.hpp"
+#include "Helpers/defs.hpp"
+#include "LinkedCell/LinkedCell_View.hpp"
 #include "LinkedCell/IPointCloud.hpp"
 #include "LinkedCell/PointCloudAdaptor.hpp"
@@ -107 +109 @@
   }
 
+  /** Determines the maximum of all element::CovalentRadius for elements present in \a &Set.
+   *
+   * I.e. the function returns a sensible cutoff criteria for bond recognition,
+   * e.g. to be used for LinkedCell_deprecated or others.
+   *
+   * \param &Set AtomSetMixin with all particles to consider
+   */
+  template <class container_type,
+            class iterator_type,
+            class const_iterator_type>
+  double getMaxPossibleBondDistance(
+      const element * const Walker,
+      const AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
+  {
+    double max_distance = 0.;
+    // get all elements
+    std::set< const element *> PresentElements;
+    for(const_iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
+      PresentElements.insert( (*AtomRunner)->getType() );
+    }
+    // create all element combinations
+    for (std::set< const element *>::const_iterator iter = PresentElements.begin();
+        iter != PresentElements.end();
+        ++iter) {
+      const range<double> MinMaxDistance(getMinMaxDistance((*iter),Walker));
+      if (MinMaxDistance.last > max_distance)
+        max_distance = MinMaxDistance.last;
+    }
+    return max_distance;
+  }
+
   /** Returns bond criterion for given pair based on a bond length matrix.
    * This calls element-version of getMinMaxDistance().
@@ -128 +161 @@
       const BondedParticle * const OtherWalker) const;
 
-  /** Creates the adjacency list for a given \a Range of iterable atoms.
-   *
-   * @param Set Range with begin and end iterator
+  /** Creates an adjacency list of the molecule.
+   * Generally, we use the CSD approach to bond recognition, that is the the distance
+   * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
+   * a threshold t = 0.4 Angstroem.
+   * To make it O(N log N) the function uses the linked-cell technique as follows:
+   * The procedure is step-wise:
+   *  -# Remove every bond in list
+   *  -# go through every atom in given \a set, check the atoms therein against all possible bond partners, add bond if true
+   *  -# correct the bond degree iteratively (single->double->triple bond)
+   *  -# finally print the bond list to \a *out if desired
+   * \param &set Container with all atoms to create adjacency for
    */
   template <class container_type,
@@ -145 +186 @@
     if (counter > 1) {
       LOG(1, "STATUS: Setting max bond distance.");
-      const double max_distance = getMaxPossibleBondDistance(Set);
+      LinkedCell::LinkedCell_View LC = getLinkedCell(getMaxPossibleBondDistance(Set));
 
       LOG(1, "STATUS: Creating LinkedCell structure for given " << counter << " atoms.");
-      PointCloudAdaptor< AtomSetMixin<container_type,iterator_type> > cloud(&Set, "SetOfAtoms");
-      LinkedCell_deprecated *LC = new LinkedCell_deprecated(cloud, max_distance);
-
-      CreateAdjacency(*LC);
-      delete (LC);
+
+      Box &domain = getDomain();
+      unsigned int CurrentTime = getTime();
+
+      unsigned int BondCount = 0;
+      // go through every atom in the set (observed cause we change its bonds)
+      for(typename AtomSetMixin<container_type,iterator_type,const_iterator_type>::iterator iter = Set.begin();
+          iter != Set.end(); ++iter) {
+        const atom * const Walker = dynamic_cast<const atom *>(*iter);
+        ASSERT(Walker != NULL,
+            "BondGraph::CreateAdjacency() - TesselPoint "
+            +(*iter)->getName()+" that was not an atom retrieved from given set");
+        LOG(2, "INFO: Current Atom is " << *Walker << ".");
+
+        // obtain all possible neighbors
+        LinkedCell::LinkedList ListOfNeighbors = LC.getAllNeighbors(
+            getMaxPossibleBondDistance(Walker->getType(), Set),
+            Walker->getPosition());
+        if (!ListOfNeighbors.empty()) {
+          // we have some possible candidates, go through each
+          for (LinkedCell::LinkedList::const_iterator neighboriter = ListOfNeighbors.begin();
+              neighboriter != ListOfNeighbors.end();
+              ++neighboriter) {
+            if ((*neighboriter) > Walker) {  // just to not add bonds from both sides
+              const atom * const OtherWalker = dynamic_cast<const atom *>(*neighboriter);
+              ASSERT(OtherWalker != NULL,
+                  "BondGraph::CreateAdjacency() - TesselPoint "
+                  +(*neighboriter)->getName()+" that was not an atom retrieved from LinkedList");
+              LOG(3, "INFO: Current other atom is " << *OtherWalker << ".");
+
+              const range<double> MinMaxDistanceSquared(
+                  getMinMaxDistanceSquared(Walker, OtherWalker));
+              const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
+              LOG(3, "INFO: Checking squared distance " << distance << " against typical bond length of " << MinMaxDistanceSquared << ".");
+              const bool status = MinMaxDistanceSquared.isInRange(distance);
+              if (status) { // create bond if distance is smaller
+                LOG(1, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
+                // directly use iter to avoid const_cast'ing Walker, too
+                //const bond * Binder =
+                const_cast<atom *>(Walker)->addBond(CurrentTime, const_cast<atom *>(OtherWalker));
+                ++BondCount;
+              } else {
+                LOG(2, "REJECT: Squared distance "
+                    << distance << " is out of squared covalent bounds "
+                    << MinMaxDistanceSquared << ".");
+              }
+
+            } else {
+              LOG(5, "REJECT: Not Adding: Wrong order.");
+            }
+          }
+        }
+      }
+      LOG(1, "I detected " << BondCount << " bonds in the molecule.");
 
       // correct bond degree by comparing valence and bond degree
@@ -242 +332 @@
     LOG(1, "STATUS: "<< bondcounter << " bonds have been parsed.");
   }
-
-  /** Creates an adjacency list of the molecule.
-   * Generally, we use the CSD approach to bond recognition, that is the the distance
-   * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
-   * a threshold t = 0.4 Angstroem.
-   * To make it O(N log N) the function uses the linked-cell technique as follows:
-   * The procedure is step-wise:
-   *  -# Remove every bond in list
-   *  -# Count the atoms in the molecule with CountAtoms()
-   *  -# partition cell into smaller linked cells of size \a bonddistance
-   *  -# put each atom into its corresponding cell
-   *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
-   *  -# correct the bond degree iteratively (single->double->triple bond)
-   *  -# finally print the bond list to \a *out if desired
-   * \param &LC Linked Cell Container with all atoms
-   */
-  void CreateAdjacency(LinkedCell_deprecated &LC) const;
 
   /** Removes all bonds within the given set of iterable atoms.
@@ -325 +398 @@
 
 private:
+  /** Convenience function to place access to World::getLinkedCell() into source module.
+   *
+   * @return ref to LinkedCell_View
+   */
+  LinkedCell::LinkedCell_View getLinkedCell(const double max_distance) const;
+
+  /** Convenience function to place access to World::getDomain() into source module.
+   *
+   * @return ref to Box
+   */
+  Box &getDomain() const;
+
+  /** Convenience function to place access to WorldTime::getTime() into source module.
+   *
+   * @return current time step
+   */
+  unsigned int getTime() const;
 
   /** Returns the BondLengthMatrix entry for a given index pair.