Changeset 1024cb for src/molecule_fragmentation.cpp

Timestamp:

May 31, 2010, 5:32:27 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, Candidate_v1.7.0, 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:

e08c46

Parents:

42af9e (diff), a7b761b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Frederik Heber <heber@…> (05/31/10 17:29:30)

git-committer:

Frederik Heber <heber@…> (05/31/10 17:32:27)

Message:

Merge commit 'jupiter/MoleculeStartEndSwitch' into CommandLineActionMapping

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/helpers.hpp
molecuilder/src/molecule.cpp
molecuilder/src/molecule_dynamics.cpp
molecuilder/src/molecule_fragmentation.cpp
molecuilder/src/molecule_geometry.cpp
molecuilder/src/molecule_graph.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/unittests/AnalysisCorrelationToPointUnitTest.cpp
molecuilder/src/unittests/listofbondsunittest.cpp

Integration of MoleculeStartEndSwitch had the following consequences:

• no more AtomCount -> getAtomCount()
• no more start/end -> begin(), end() and iterator
• no more decent ordering in atomic ids (hence, Simple_configuration/8 and Domain/5, Domain/6 now check by comparing sorted xyz, not confs)

There is still a huge problem with bonds. One test runs into an endless loop.

Signed-off-by: Frederik Heber <heber@…>

File:

• src/molecule_fragmentation.cpp (modified) (25 diffs)

src/molecule_fragmentation.cpp

@@ -39 +39 @@
   int FragmentCount;
   // get maximum bond degree
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    c = (Walker->ListOfBonds.size() > c) ? Walker->ListOfBonds.size() : c;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c;
   }
   FragmentCount = NoNonHydrogen*(1 << (c*order));
@@ -353 +351 @@
 map<double, pair<int,int> >  * ReMapAdaptiveCriteriaListToValue(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
 {
-  atom *Walker = mol->start;
+  atom *Walker = NULL;
   map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
   DoLog(1) && (Log() << Verbose(1) << "Root candidate list is: " << endl);
@@ -385 +383 @@
 bool MarkUpdateCandidates(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
 {
-  atom *Walker = mol->start;
+  atom *Walker = NULL;
   int No = -1;
   bool status = false;
@@ -429 +427 @@
 bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
 {
-  atom *Walker = start;
   bool status = false;
 
   // initialize mask list
-  for(int i=AtomCount;i--;)
+  for(int i=getAtomCount();i--;)
     AtomMask[i] = false;
 
   if (Order < 0) { // adaptive increase of BondOrder per site
-    if (AtomMask[AtomCount] == true)  // break after one step
+    if (AtomMask[getAtomCount()] == true)  // break after one step
       return false;
 
@@ -451 +448 @@
     if (AdaptiveCriteriaList->empty()) {
       DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
-      while (Walker->next != end) {
-        Walker = Walker->next;
+      for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     #ifdef ADDHYDROGEN
-        if (Walker->type->Z != 1) // skip hydrogen
+        if ((*iter)->type->Z != 1) // skip hydrogen
     #endif
         {
-          AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
+          AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
           status = true;
         }
@@ -472 +468 @@
     delete[](FinalRootCandidates);
   } else { // global increase of Bond Order
-    while (Walker->next != end) {
-      Walker = Walker->next;
+    for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
   #ifdef ADDHYDROGEN
-      if (Walker->type->Z != 1) // skip hydrogen
+      if ((*iter)->type->Z != 1) // skip hydrogen
   #endif
       {
-        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
-        if ((Order != 0) && (Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
+        AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
+        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->nr]))
           status = true;
       }
     }
-    if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
+    if ((!Order) && (!AtomMask[getAtomCount()]))  // single stepping, just check
       status = true;
 
@@ -494 +489 @@
   }
 
-  PrintAtomMask(AtomMask, AtomCount); // for debugging
+  PrintAtomMask(AtomMask, getAtomCount()); // for debugging
 
   return status;
@@ -510 +505 @@
     return false;
   }
-  SortIndex = new int[AtomCount];
-  for(int i=AtomCount;i--;)
+  SortIndex = new int[getAtomCount()];
+  for(int i=getAtomCount();i--;)
     SortIndex[i] = -1;
 
@@ -518 +513 @@
 
   return true;
+};
+
+
+
+/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
+ * \param *start begin of list (STL iterator, i.e. first item)
+ * \paran *end end of list (STL iterator, i.e. one past last item)
+ * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
+ * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
+ * \return true - success, false - failure
+ */
+bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
+{
+  bool status = true;
+  int AtomNo;
+
+  if (LookupTable != NULL) {
+    Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
+    return false;
+  }
+
+  // count them
+  if (count == 0) {
+    for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
+      count = (count < (*iter)->GetTrueFather()->nr) ? (*iter)->GetTrueFather()->nr : count;
+    }
+  }
+  if (count <= 0) {
+    Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
+    return false;
+  }
+
+  // allocate and fill
+  LookupTable = new atom *[count];
+  if (LookupTable == NULL) {
+    eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
+    performCriticalExit();
+    status = false;
+  } else {
+    for (int i=0;i<count;i++)
+      LookupTable[i] = NULL;
+    for (molecule::iterator iter = begin(); iter != end(); ++iter) {
+      AtomNo = (*iter)->GetTrueFather()->nr;
+      if ((AtomNo >= 0) && (AtomNo < count)) {
+        //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl;
+        LookupTable[AtomNo] = (*iter);
+      } else {
+        Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl;
+        status = false;
+        break;
+      }
+    }
+  }
+
+  return status;
 };
 
@@ -541 +591 @@
 {
   MoleculeListClass *BondFragments = NULL;
-  int *MinimumRingSize = new int[AtomCount];
+  int *MinimumRingSize = new int[getAtomCount()];
   int FragmentCounter;
   MoleculeLeafClass *MolecularWalker = NULL;
@@ -569 +619 @@
 
   // create lookup table for Atom::nr
-  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(start, end, ListOfAtoms, AtomCount);
+  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, getAtomCount());
 
   // === compare it with adjacency file ===
@@ -579 +629 @@
 
   // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
-  for(int i=AtomCount;i--;)
-    MinimumRingSize[i] = AtomCount;
+  for(int i=getAtomCount();i--;)
+    MinimumRingSize[i] = getAtomCount();
   MolecularWalker = Subgraphs;
   FragmentCounter = 0;
@@ -586 +636 @@
     MolecularWalker = MolecularWalker->next;
     // fill the bond structure of the individually stored subgraphs
-  MolecularWalker->FillBondStructureFromReference(this, FragmentCounter, ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
+    MolecularWalker->FillBondStructureFromReference(this, FragmentCounter, ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
     DoLog(0) && (Log() << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
     LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
 //    // check the list of local atoms for debugging
 //    Log() << Verbose(0) << "ListOfLocalAtoms for this subgraph is:" << endl;
-//    for (int i=0;i<AtomCount;i++)
+//    for (int i=0;i<getAtomCount();i++)
 //      if (ListOfLocalAtoms[FragmentCounter][i] == NULL)
 //        Log() << Verbose(0) << "\tNULL";
@@ -617 +667 @@
   // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
   KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
-  AtomMask = new bool[AtomCount+1];
-  AtomMask[AtomCount] = false;
+  AtomMask = new bool[getAtomCount()+1];
+  AtomMask[getAtomCount()] = false;
   FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
   while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
     FragmentationToDo = FragmentationToDo || CheckOrder;
-    AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
+    AtomMask[getAtomCount()] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
     // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
     Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0));
@@ -762 +812 @@
 bool molecule::ParseOrderAtSiteFromFile(char *path)
 {
-  unsigned char *OrderArray = new unsigned char[AtomCount];
-  bool *MaxArray = new bool[AtomCount];
+  unsigned char *OrderArray = new unsigned char[getAtomCount()];
+  bool *MaxArray = new bool[getAtomCount()];
   bool status;
   int AtomNr, value;
@@ -769 +819 @@
   ifstream file;
 
-  for(int i=0;i<AtomCount;i++) {
+  for(int i=0;i<getAtomCount();i++) {
     OrderArray[i] = 0;
     MaxArray[i] = false;
@@ -871 +921 @@
   atom *OtherFather = NULL;
   atom *FatherOfRunner = NULL;
-  Leaf->CountAtoms();
-
-  atom *Runner = Leaf->start;
-  while (Runner->next != Leaf->end) {
-    Runner = Runner->next;
+
+#ifdef ADDHYDROGEN
+  molecule::const_iterator runner;
+#endif
+  // we increment the iter just before skipping the hydrogen
+  for (molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end();) {
     LonelyFlag = true;
-    FatherOfRunner = Runner->father;
+    FatherOfRunner = (*iter)->father;
+    ASSERT(FatherOfRunner,"Atom without father found");
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
@@ -888 +940 @@
 //            Log() << Verbose(3) << "Adding Bond: ";
 //            Log() << Verbose(0) <<
-            Leaf->AddBond(Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree);
+            Leaf->AddBond((*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree);
 //            Log() << Verbose(0) << "." << endl;
-            //NumBonds[Runner->nr]++;
+            //NumBonds[(*iter)->nr]++;
           } else {
 //            Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
@@ -898 +950 @@
 //          Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl;
 #ifdef ADDHYDROGEN
-          //Log() << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), Runner, FatherOfRunner, OtherFather, IsAngstroem))
+          //Log() << Verbose(3) << "Adding Hydrogen to " << (*iter)->Name << " and a bond in between." << endl;
+          if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
             exit(1);
 #endif
-          //NumBonds[Runner->nr] += Binder->BondDegree;
+          //NumBonds[(*iter)->nr] += Binder->BondDegree;
         }
       }
     } else {
-      DoeLog(1) && (eLog()<< Verbose(1) << "Son " << Runner->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
-    }
-    if ((LonelyFlag) && (Leaf->AtomCount > 1)) {
-      DoLog(0) && (Log() << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl);
-    }
+    DoeLog(1) && (eLog()<< Verbose(1) << "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
+    }
+    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
+      DoLog(0) && (Log() << Verbose(0) << **iter << "has got bonds only to hydrogens!" << endl);
+    }
+    ++iter;
 #ifdef ADDHYDROGEN
-    while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
-      Runner = Runner->next;
+    while ((iter != Leaf->end()) && ((*iter)->type->Z == 1)){ // skip added hydrogen
+      iter++;
+    }
 #endif
   }
@@ -928 +982 @@
 molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
 {
-  atom **SonList = new atom*[AtomCount];
+  atom **SonList = new atom*[getAtomCount()];
   molecule *Leaf = World::getInstance().createMolecule();
 
-  for(int i=0;i<AtomCount;i++)
+  for(int i=0;i<getAtomCount();i++)
     SonList[i] = NULL;
 
@@ -1561 +1615 @@
   FragmentSearch.FragmentSet = new KeySet;
   FragmentSearch.Root = FindAtom(RootKeyNr);
-  FragmentSearch.ShortestPathList = new int[AtomCount];
-  for (int i=AtomCount;i--;) {
+  FragmentSearch.ShortestPathList = new int[getAtomCount()];
+  for (int i=getAtomCount();i--;) {
     FragmentSearch.ShortestPathList[i] = -1;
   }
 
   // Construct the complete KeySet which we need for topmost level only (but for all Roots)
-  atom *Walker = start;
   KeySet CompleteMolecule;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    CompleteMolecule.insert(Walker->GetTrueFather()->nr);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    CompleteMolecule.insert((*iter)->GetTrueFather()->nr);
   }
 
@@ -1582 +1634 @@
     RootKeyNr = RootStack.front();
     RootStack.pop_front();
-    Walker = FindAtom(RootKeyNr);
+    atom *Walker = FindAtom(RootKeyNr);
     // check cyclic lengths
     //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
@@ -1671 +1723 @@
   Vector Translationvector;
   //class StackClass<atom *> *CompStack = NULL;
-  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
+  class StackClass<atom *> *AtomStack = new StackClass<atom *>(getAtomCount());
   bool flag = true;
 
   DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl);
 
-  ColorList = new enum Shading[AtomCount];
-  for (int i=0;i<AtomCount;i++)
+  ColorList = new enum Shading[getAtomCount()];
+  for (int i=0;i<getAtomCount();i++)
     ColorList[i] = (enum Shading)0;
   while (flag) {
     // remove bonds that are beyond bonddistance
-    for(int i=NDIM;i--;)
-      Translationvector[i] = 0.;
+    Translationvector.Zero();
     // scan all bonds
     Binder = first;
@@ -1711 +1762 @@
       Log() << Verbose(0) << Translationvector <<  endl;
       // apply to all atoms of first component via BFS
-      for (int i=AtomCount;i--;)
+      for (int i=getAtomCount();i--;)
         ColorList[i] = white;
       AtomStack->Push(Binder->leftatom);
@@ -1735 +1786 @@
     //delete(CompStack);
   }
-
   // free allocated space from ReturnFullMatrixforSymmetric()
   delete(AtomStack);