Changeset 41baaf for src

Timestamp:

Jul 30, 2008, 1:11:30 PM (16 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:

af6d8f

Parents:

683914

Message:

2 BUGFIXES: CreateAdjacencyList() and CyclicStructureAnalysis()

CreateAdjacencyList(): with benzophenone and indigo it was found that the correction of the bond degree was not working properly. Carbon bound to oxygen and other carbons would not raise the bond degree to oxygen preferentially over the carbon bonds, which lead to wrong bond degrees in the end. We refrained from using matrix schemes as these are not O(N) in a simple manner, but remained with the local scheme that now checks each also mismatching bond partner for its number of bonds and takes preferentially the one with the fewest bonds as the increase candidate.
CyclicStructureAnalysis(): with biphenylene and cholesterol it was found that the cycle detection did not work as expected. Cycles were found twice due to a back edge being shared by two cycles at the same time (and it's randomly dependend on the vertex ids which cycle is found first). Hence, we now do a backtracking of a cycle candidate, looking whether each of its vertices is not already set to be cyclic. If so, we removed the last added vertices from lists and stacks (i.e. the closing edges), so that the other BFS finger with the real candidate may still reach the looked-for backedge (otherwise the vertex would be marked gray and hence not visited).

File:

• src/molecules.cpp (modified) (7 diffs)

src/molecules.cpp

@@ -1442 +1442 @@
 void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
 {
-  atom *Walker = NULL, *OtherWalker = NULL;
-  int No, NoBonds;
+  atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
+  int No, NoBonds, CandidateBondNo;
   int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
   molecule **CellList;
@@ -1556 +1556 @@
     CreateListOfBondsPerAtom(out);
                 
-    // correct Bond degree of each bond by checking of updated(!) sum of bond degree for an atom match its valence count
-    // bond degrres are correctled iteratively by one, so that 2-2 instead of 1-3 or 3-1 corrections are favoured: We want
-    // a rather symmetric distribution of higher bond degrees
+    // 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;
@@ -1567 +1568 @@
         while (Walker->next != end) { // go through every atom
           Walker = Walker->next;
-          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 first partner (updated!), might have changed during last bond partner
-                  NoBonds = 0;
-            for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
-              NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
-                  *out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-                  if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check NoBonds of other atom
+          // 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;
+            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: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-                    if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) // increase bond degree by one
-                ListOfBondsPerAtom[Walker->nr][i]->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) {
+              ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
+              *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
+            }
           }
                     }
@@ -1801 +1812 @@
 
 /** Analyses the cycles found and returns minimum of all cycle lengths.
+ * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root,
+ * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
+ * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
+ * as cyclic and print out the cycles. 
  * \param *out output stream for debugging
  * \param *BackEdgeStack stack with all back edges found during DFS scan
@@ -1812 +1827 @@
   enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
   class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring 
-  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms 
+  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
   bond *Binder = NULL, *BackEdge = NULL;
@@ -1843 +1858 @@
     BFSStack->Push(Walker);
     TouchedStack->Push(Walker);
-    //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
+    *out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
     OtherAtom = NULL;
-    while ((Walker != Root) && ((OtherAtom == NULL) || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]))) {  // look for Root
+    do {  // look for Root
       Walker = BFSStack->PopFirst();
-      //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
+      *out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
       for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
         Binder = ListOfBondsPerAtom[Walker->nr][i];
         if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
           OtherAtom = Binder->GetOtherAtom(Walker);
-          //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
-          if (ColorList[OtherAtom->nr] == white) {
-            TouchedStack->Push(OtherAtom);
-            ColorList[OtherAtom->nr] = lightgray;
-            PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
-            ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-            //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
-            if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
-              //*out << Verbose(3) << "Putting OtherAtom into queue." << endl;
-              BFSStack->Push(OtherAtom);
+#ifdef ADDHYDROGEN          
+          if (OtherAtom->type->Z != 1) {
+#endif
+            *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
+            if (ColorList[OtherAtom->nr] == white) {
+              TouchedStack->Push(OtherAtom);
+              ColorList[OtherAtom->nr] = lightgray;
+              PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
+              ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
+              *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+              //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
+                *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
+                BFSStack->Push(OtherAtom);
+              //}
+            } else {
+              *out << Verbose(3) << "Not Adding, has already been visited." << endl;
             }
+            if (OtherAtom == Root)
+              break;
+#ifdef ADDHYDROGEN          
           } else {
-            //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
+            *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
+            ColorList[OtherAtom->nr] = black;
           }
+#endif
         } else {
-          //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
+          *out << Verbose(2) << "Bond " << *Binder << " not Visiting, is the back edge." << endl;
         }
       }
       ColorList[Walker->nr] = black;
-      //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
-    }
+      *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
+      if (OtherAtom == Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
+        // step through predecessor list
+        while (OtherAtom != BackEdge->rightatom) {
+          if (!OtherAtom->GetTrueFather()->IsCyclic)  // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
+            break;
+          else
+            OtherAtom = PredecessorList[OtherAtom->nr];
+        }
+        if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
+          *out << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl;\
+          do {
+            OtherAtom = TouchedStack->PopLast();
+            if (PredecessorList[OtherAtom->nr] == Walker) {
+              *out << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl;
+              PredecessorList[OtherAtom->nr] = NULL;
+              ShortestPathList[OtherAtom->nr] = -1;
+              ColorList[OtherAtom->nr] = white;
+              BFSStack->RemoveItem(OtherAtom);
+            }
+          } while ((!TouchedStack->IsEmpty()) && (PredecessorList[OtherAtom->nr] == NULL));
+          TouchedStack->Push(OtherAtom);  // last was wrongly popped
+          OtherAtom = BackEdge->rightatom; // set to not Root
+        } else
+          OtherAtom = Root;
+      }
+    } while ((OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
     
-    if (Walker == Root) {
+    if (OtherAtom == Root) {
       // now climb back the predecessor list and thus find the cycle members
       NumCycles++;
@@ -1880 +1931 @@
       Root->GetTrueFather()->IsCyclic = true;
       *out << Verbose(1) << "Found ring contains: ";
+      Walker = Root;
       while (Walker != BackEdge->rightatom) {
         *out << Walker->Name << " <-> ";