Changeset b8b75d

Timestamp:

Oct 17, 2009, 7:47:58 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:

266237

Parents:

07051c

Message:

Refactored CreateAdjacencyList().

• new auxiliary functions: molecule::CorrectBondDegree(), molecule::OutputBondsList() and molecule::CountAtomsBonds()
• CreateAdjacencyList() makes use of the above and the bond finding is implemented now via LinkedCell list, however with an additional TesselPoint::nr to atom* list (which is just a temporary fix until a thought comes up how to have LinkedCell of atom * in the first place).

Location:

src

Files:

• molecule.hpp (modified) (1 diff)
• molecule_graph.cpp (modified) (6 diffs)

src/molecule.hpp

@@ -248 +248 @@
   void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
   void CreateListOfBondsPerAtom(ofstream *out);
+  int CorrectBondDegree(ofstream *out);
+  void OutputBondsList(ofstream *out);
+  int CountAtomsBonds(int nr);
+
 
   // Graph analysis

src/molecule_graph.cpp

@@ -11 +11 @@
 #include "element.hpp"
 #include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "lists.hpp"
 #include "memoryallocator.hpp"
@@ -73 +74 @@
 {
 
-  atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
-  int No, NoBonds, CandidateBondNo;
-  int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
-  molecule **CellList;
+  atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  atom **AtomMap = NULL;
+  int n[NDIM];
   double distance, MinDistance, MaxDistance;
-  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
-  Vector x;
-  int FalseBondDegree = 0;
+  LinkedCell *LC = NULL;
+  LinkedNodes *List = NULL;
+  LinkedNodes *OtherList = NULL;
 
   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
@@ -94 +95 @@
 
   if (AtomCount != 0) {
-    // 1. find divisor for each axis, such that a sphere with radius of at least bonddistance can be placed into each cell
-    j=-1;
-    for (int i=0;i<NDIM;i++) {
-      j += i+1;
-      divisor[i] = (int)floor(cell_size[j]/bonddistance); // take smaller value such that size of linked cell is at least bonddistance
-      //*out << Verbose(1) << "divisor[" << i << "]  = " << divisor[i] << "." << endl;
-    }
-    // 2a. allocate memory for the cell list
-    NumberCells = divisor[0]*divisor[1]*divisor[2];
-    *out << Verbose(1) << "Allocating " << NumberCells << " cells." << endl;
-    CellList = Malloc<molecule*>(NumberCells, "molecule::CreateAdjacencyList - ** CellList");
-    for (int i=NumberCells;i--;)
-      CellList[i] = NULL;
-
-    // 2b. put all atoms into its corresponding list
+    LC = new LinkedCell(this, bonddistance);
+
+    // create a list to map Tesselpoint::nr to atom *
+    AtomMap = Malloc<atom *>(AtomCount, "molecule::CreateAdjacencyList - **AtomCount");
     Walker = start;
-    while(Walker->next != end) {
+    while (Walker->next != end) {
       Walker = Walker->next;
-      //*out << Verbose(1) << "Current atom is " << *Walker << " with coordinates ";
-      //Walker->x.Output(out);
-      //*out << "." << endl;
-      // compute the cell by the atom's coordinates
-      j=-1;
-      for (int i=0;i<NDIM;i++) {
-        j += i+1;
-        x.CopyVector(&(Walker->x));
-        x.KeepPeriodic(out, matrix);
-        n[i] = (int)floor(x.x[i]/cell_size[j]*(double)divisor[i]);
-      }
-      index = n[2] + (n[1] + n[0] * divisor[1]) * divisor[2];
-      //*out << Verbose(1) << "Atom " << *Walker << " goes into cell number [" << n[0] << "," << n[1] << "," << n[2] << "] = " << index << "." << endl;
-      // add copy atom to this cell
-      if (CellList[index] == NULL)  // allocate molecule if not done
-        CellList[index] = new molecule(elemente);
-      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
-      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
-    }
-    //for (int i=0;i<NumberCells;i++)
-      //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
-
+      AtomMap[Walker->nr] = Walker;
+    }
 
     // 3a. go through every cell
-    for (N[0]=divisor[0];N[0]--;)
-      for (N[1]=divisor[1];N[1]--;)
-        for (N[2]=divisor[2];N[2]--;) {
-          Index = N[2] + (N[1] + N[0] * divisor[1]) * divisor[2];
-          if (CellList[Index] != NULL) { // if there atoms in this cell
-            //*out << Verbose(1) << "Current cell is " << Index << "." << endl;
-            // 3b. for every atom therein
-            Walker = CellList[Index]->start;
-            while (Walker->next != CellList[Index]->end) {  // go through every atom
-              Walker = Walker->next;
+    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]++) {
+          List = LC->GetCurrentCell();
+          //*out << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+          if (List != NULL) {
+            for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+              Walker = AtomMap[(*Runner)->nr];
               //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
               // 3c. check for possible bond between each atom in this and every one in the 27 cells
@@ -151 +119 @@
                 for (n[1]=-1;n[1]<=1;n[1]++)
                   for (n[2]=-1;n[2]<=1;n[2]++) {
-                     // compute the index of this comparison cell and make it periodic
-                    index = ((N[2]+n[2]+divisor[2])%divisor[2]) + (((N[1]+n[1]+divisor[1])%divisor[1]) + ((N[0]+n[0]+divisor[0])%divisor[0]) * divisor[1]) * divisor[2];
-                    //*out << Verbose(1) << "Number of comparison cell is " << index << "." << endl;
-                    if (CellList[index] != NULL) {  // if there are any atoms in this cell
-                      OtherWalker = CellList[index]->start;
-                      while(OtherWalker->next != CellList[index]->end) {  // go through every atom in this cell
-                        OtherWalker = OtherWalker->next;
-                        //*out << Verbose(0) << "Current comparison atom is " << *OtherWalker << "." << endl;
-                        /// \todo periodic check is missing here!
-                        //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
-                        MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
-                        MinDistance *= (IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem;
-                        MaxDistance = MinDistance + BONDTHRESHOLD;
-                        MinDistance -= BONDTHRESHOLD;
-                        distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
-                        if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
-                          //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
-                          AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
-                        } else {
-                          //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
+                    OtherList = LC->GetRelativeToCurrentCell(n);
+                    //*out << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+                    if (OtherList != NULL) {
+                      for (LinkedNodes::iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
+                        if ((*OtherRunner)->nr > Walker->nr) {
+                          OtherWalker = AtomMap[(*OtherRunner)->nr];
+                          //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
+                          MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
+                          MinDistance *= (IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem;
+                          MaxDistance = MinDistance + BONDTHRESHOLD;
+                          MinDistance -= BONDTHRESHOLD;
+                          distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
+                          if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
+                            //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
+                            AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
+                          } else {
+                            //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
+                          }
                         }
                       }
@@ -178 +144 @@
           }
         }
-
-
-
-    // 4. free the cell again
-    for (int i=NumberCells;i--;)
-      if (CellList[i] != NULL) {
-        delete(CellList[i]);
-      }
-    Free(&CellList);
+    Free(&AtomMap);
+    delete(LC);
+    *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl;
 
     // create the adjacency list per atom
     CreateListOfBondsPerAtom(out);
 
-    // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
-    // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
-    // preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
-    // double bonds as was expected.
-    if (BondCount != 0) {
-      NoCyclicBonds = 0;
-      *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
-      do {
-        No = 0; // No acts as breakup flag (if 1 we still continue)
-        Walker = start;
-        while (Walker->next != end) { // go through every atom
-          Walker = Walker->next;
-          // count valence of first partner
-          NoBonds = 0;
-          for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
-            NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-          *out << Verbose(3) << "Walker " << *Walker << ": " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-          if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
-            Candidate = NULL;
-            CandidateBondNo = -1;
-            for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
-              OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-              // count valence of second partner
-              NoBonds = 0;
-              for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
-                NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
-              *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-              if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
-                if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
-                  Candidate = OtherWalker;
-                  CandidateBondNo = i;
-                  *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
-                }
-              }
-            }
-            if ((Candidate != NULL) && (CandidateBondNo != -1)) {
-              ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
-              *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
-            } else
-              *out << Verbose(2) << "Could not find correct degree for atom " << *Walker << "." << endl;
-              FalseBondDegree++;
-          }
-        }
-      } while (No);
-    *out << " done." << endl;
-    } else
-      *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
-    *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
+    // correct bond degree by comparing valence and bond degree
+    CorrectBondDegree(out);
 
     // output bonds for debugging (if bond chain list was correctly installed)
-    *out << Verbose(1) << endl << "From contents of bond chain list:";
-    bond *Binder = first;
-    while(Binder->next != last) {
-      Binder = Binder->next;
-      *out << *Binder << "\t" << endl;
-    }
-    *out << endl;
+    OutputBondsList(out);
   } else
     *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
-  Free(&matrix);
-
+};
+
+/** Prints a list of all bonds to \a *out.
+ * \param output stream
+ */
+void molecule::OutputBondsList(ofstream *out)
+{
+  *out << Verbose(1) << endl << "From contents of bond chain list:";
+  bond *Binder = first;
+  while(Binder->next != last) {
+    Binder = Binder->next;
+    *out << *Binder << "\t" << endl;
+  }
+  *out << endl;
+};
+
+/** Sums up the number of bonds times bond::BondDegree the atom with index \a nr has.
+ * \param nr indexof atom
+ * \return number of bonds each weighted with its bond::BondDegree
+ */
+int molecule::CountAtomsBonds(int nr)
+{
+  int NoBonds = 0;
+  for(int j=0;j<NumberOfBondsPerAtom[nr];j++)
+    NoBonds += ListOfBondsPerAtom[nr][j]->BondDegree;
+  return NoBonds;
+};
+
+/** correct bond degree by comparing valence and bond degree.
+ * correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
+ * iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
+ * preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
+ * double bonds as was expected.
+ * \param *out output stream for debugging
+ * \return number of bonds that could not be corrected
+ */
+int molecule::CorrectBondDegree(ofstream *out)
+{
+  int No = 0;
+  int NoBonds = 0;
+  int CandidateBondNo = 0;
+  int FalseBondDegree = 0;
+  atom *Walker = NULL;
+  atom *Candidate = NULL;
+  atom *OtherWalker = NULL;
+
+  if (BondCount != 0) {
+    NoCyclicBonds = 0;
+    *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
+    do {
+      No = 0; // No acts as breakup flag (if 1 we still continue)
+      Walker = start;
+      while (Walker->next != end) { // go through every atom
+        Walker = Walker->next;
+
+        NoBonds = CountAtomsBonds(Walker->nr);
+        *out << Verbose(3) << "Walker " << *Walker << ": " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+        if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
+          Candidate = NULL;
+          CandidateBondNo = -1;
+          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
+            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+            NoBonds = CountAtomsBonds(OtherWalker->nr);
+            *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+            if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
+              if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
+                Candidate = OtherWalker;
+                CandidateBondNo = i;
+                *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
+              }
+            }
+          }
+          if ((Candidate != NULL) && (CandidateBondNo != -1)) {
+            ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
+            *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
+          } else {
+            *out << Verbose(2) << "Could not find correct degree for atom " << *Walker << "." << endl;
+            FalseBondDegree++;
+          }
+        }
+      }
+    } while (No);
+    *out << " done." << endl;
+  } else {
+    *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
+  }
+  *out << FalseBondDegree << " bonds could not be corrected." << endl;
+
+  return (FalseBondDegree);
 };
 
@@ -285 +280 @@
   return No;
 };
+
+
 /** Returns Shading as a char string.
  * \param color the Shading