Changes in src/moleculelist.cpp [f66195:c9bce3e]

File:

• src/moleculelist.cpp (modified) (18 diffs)

src/moleculelist.cpp

@@ -12 +12 @@
 #include "helpers.hpp"
 #include "linkedcell.hpp"
+#include "lists.hpp"
 #include "molecule.hpp"
 #include "memoryallocator.hpp"
@@ -306 +307 @@
 bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
 {
+  LinkedCell *LCList = NULL;
+  Tesselation *TesselStruct = NULL;
   if ((srcmol == NULL) || (mol == NULL)) {
     cout << Verbose(1) << "ERROR: Either fixed or variable molecule is given as NULL." << endl;
@@ -312 +315 @@
 
   // calculate envelope for *mol
-  LinkedCell *LCList = new LinkedCell(mol, 8.);
-  FindNonConvexBorder((ofstream *)&cout, mol, LCList, 4., NULL);
-  if (mol->TesselStruct == NULL) {
+  LCList = new LinkedCell(mol, 8.);
+  FindNonConvexBorder((ofstream *)&cout, mol, TesselStruct, (const LinkedCell *&)LCList, 4., NULL);
+  if (TesselStruct == NULL) {
     cout << Verbose(1) << "ERROR: Could not tesselate the fixed molecule." << endl;
     return false;
   }
   delete(LCList);
-  LCList = new LinkedCell(mol->TesselStruct, 8.);  // re-create with boundary points only!
+  LCList = new LinkedCell(TesselStruct, 8.);  // re-create with boundary points only!
 
   // prepare index list for bonds
@@ -333 +336 @@
     Walker = Walker->next;
     cout << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl;
-    if (!mol->TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
+    if (!TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
       CopyAtoms[Walker->nr] = new atom(Walker);
       mol->AddAtom(CopyAtoms[Walker->nr]);
@@ -376 +379 @@
   atom *Walker = NULL;
   atom *Runner = NULL;
+  bond *Binder = NULL;
   double ***FitConstant = NULL, **correction = NULL;
   int a, b;
@@ -419 +423 @@
 
   // 0b. allocate memory for constants
-  FitConstant = Malloc<double**>(3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  FitConstant = Calloc<double**>(3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
   for (int k = 0; k < 3; k++) {
-    FitConstant[k] = Malloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    FitConstant[k] = Calloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
     for (int i = a; i--;) {
-      FitConstant[k][i] = Malloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+      FitConstant[k][i] = Calloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     }
   }
@@ -468 +472 @@
 
   // 0d. allocate final correction matrix
-  correction = Malloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **correction");
+  correction = Calloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **correction");
   for (int i = a; i--;)
-    correction[i] = Malloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
+    correction[i] = Calloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
 
   // 1a. go through every molecule in the list
@@ -482 +486 @@
     while (Walker->next != (*ListRunner)->end) {
       Walker = Walker->next;
-      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
       if ((Walker->type->Z == 1) && ((Walker->father == NULL)
           || (Walker->father->type->Z != 1))) { // if it's a hydrogen
@@ -488 +492 @@
         while (Runner->next != (*ListRunner)->end) {
           Runner = Runner->next;
-          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
-          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
+          Binder = *(Runner->ListOfBonds.begin());
+          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (Binder->GetOtherAtom(Runner) != Binder->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
             // 4. evaluate the morse potential for each matrix component and add up
             distance = Runner->x.Distance(&Walker->x);
@@ -544 +549 @@
   output.close();
   // 6. free memory of parsed matrices
-  FitConstant = Malloc<double**>(a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
   for (int k = 0; k < 3; k++) {
-    FitConstant[k] = Malloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
     for (int i = a; i--;) {
-      FitConstant[k][i] = Malloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
-    }
-  }
+      Free(&FitConstant[k][i]);
+    }
+    Free(&FitConstant[k]);
+  }
+  Free(&FitConstant);
   cout << "done." << endl;
   return true;
@@ -708 +713 @@
     //outputFragment.close();
     //outputFragment.clear();
-    delete (FragmentNumber);
-    //Free(&FragmentNumber);
+    Free(&FragmentNumber);
   }
   cout << " done." << endl;
@@ -730 +734 @@
 };
 
+/** Dissects given \a *mol into connected subgraphs and inserts them as new molecules but with old atoms into \a this.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms to dissect
+ * \param *configuration config with BondGraph
+ */
+void MoleculeListClass::DissectMoleculeIntoConnectedSubgraphs(ofstream * const out, molecule * const mol, config * const configuration)
+{
+  // 1. dissect the molecule into connected subgraphs
+  configuration ->BG->ConstructBondGraph(out, mol);
+
+  // 2. scan for connected subgraphs
+  MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
+  class StackClass<bond *> *BackEdgeStack = NULL;
+  Subgraphs = mol->DepthFirstSearchAnalysis(out, BackEdgeStack);
+  delete(BackEdgeStack);
+
+  // 3. dissect (the following construct is needed to have the atoms not in the order of the DFS, but in
+  // the original one as parsed in)
+  // TODO: Optimize this, when molecules just contain pointer list of global atoms!
+
+  // 4a. create array of molecules to fill
+  const int MolCount = Subgraphs->next->Count();
+  molecule **molecules = Malloc<molecule *>(MolCount, "config::Load() - **molecules");
+  for (int i=0;i<MolCount;i++) {
+    molecules[i] = (molecule*) new molecule(mol->elemente);
+    molecules[i]->ActiveFlag = true;
+    insert(molecules[i]);
+  }
+
+  // 4b. create and fill map of which atom is associated to which connected molecule (note, counting starts at 1)
+  int FragmentCounter = 0;
+  int *MolMap = Calloc<int>(mol->AtomCount, "config::Load() - *MolMap");
+  MoleculeLeafClass *MolecularWalker = Subgraphs;
+  atom *Walker = NULL;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    Walker = MolecularWalker->Leaf->start;
+    while (Walker->next != MolecularWalker->Leaf->end) {
+      Walker = Walker->next;
+      MolMap[Walker->GetTrueFather()->nr] = FragmentCounter+1;
+    }
+    FragmentCounter++;
+  }
+
+  // 4c. relocate atoms to new molecules and remove from Leafs
+  Walker = mol->start;
+  while (mol->start->next != mol->end) {
+    Walker = mol->start->next;
+    if ((Walker->nr <0) || (Walker->nr >= mol->AtomCount)) {
+      cerr << "Index of atom " << *Walker << " is invalid!" << endl;
+      performCriticalExit();
+    }
+    FragmentCounter = MolMap[Walker->nr];
+    if (FragmentCounter != 0) {
+      cout << Verbose(3) << "Re-linking " << *Walker << "..." << endl;
+      unlink(Walker);
+      molecules[FragmentCounter-1]->AddAtom(Walker);    // counting starts at 1
+    } else {
+      cerr << "Atom " << *Walker << " not associated to molecule!" << endl;
+      performCriticalExit();
+    }
+  }
+  // 4d. we don't need to redo bonds, as they are connected subgraphs and still maintained their ListOfBonds
+  // 4e. free Leafs
+  MolecularWalker = Subgraphs;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    delete(MolecularWalker->previous);
+  }
+  delete(MolecularWalker);
+  Free(&MolMap);
+  Free(&molecules);
+  cout << Verbose(1) << "I scanned " << FragmentCounter << " molecules." << endl;
+};
 
 /******************************************* Class MoleculeLeafClass ************************************************/
@@ -810 +888 @@
  * \return true - success, false - faoilure
  */
-bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
-{
-  atom *Walker = NULL, *OtherWalker = NULL;
-  bond *Binder = NULL;
+bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, const molecule * const reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
+{
+  atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  atom *Father = NULL;
   bool status = true;
   int AtomNo;
@@ -825 +904 @@
 
   if (status) {
-    *out << Verbose(1) << "Creating adjacency list for subgraph " << this
-        << "." << endl;
+    *out << Verbose(1) << "Creating adjacency list for subgraph " << Leaf << "." << endl;
+    // remove every bond from the list
+    bond *Binder = NULL;
+    while (Leaf->last->previous != Leaf->first) {
+      Binder = Leaf->last->previous;
+      Binder->leftatom->UnregisterBond(Binder);
+      Binder->rightatom->UnregisterBond(Binder);
+      removewithoutcheck(Binder);
+    }
+
     Walker = Leaf->start;
     while (Walker->next != Leaf->end) {
       Walker = Walker->next;
-      AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker
-      for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all
-        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
-        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
+      Father = Walker->GetTrueFather();
+      AtomNo = Father->nr; // global id of the current walker
+      for (BondList::const_iterator Runner = Father->ListOfBonds.begin(); Runner != Father->ListOfBonds.end(); (++Runner)) {
+        OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
         if (OtherWalker != NULL) {
           if (OtherWalker->nr > Walker->nr)
-            Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
+            Leaf->AddBond(Walker, OtherWalker, (*Runner)->BondDegree);
         } else {
-          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
+          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
           status = false;
         }
       }
     }
-    Leaf->CreateListOfBondsPerAtom(out);
-    FragmentCounter++;
-    if (next != NULL)
-      status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
-    FragmentCounter--;
   }
 
@@ -856 +938 @@
       Free(&ListOfLocalAtoms);
   }
-  FragmentCounter--;
   *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
   return status;
@@ -903 +984 @@
 
 /** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
- * \param *out output stream fro debugging
+ * \param *out output stream from debugging
  * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
  * \param FragmentCounter counts the fragments as we move along the list
  * \param GlobalAtomCount number of atoms in the complete molecule
  * \param &FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
- * \return true - succes, false - failure
+ * \return true - success, false - failure
  */
 bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList)
@@ -914 +995 @@
   bool status = true;
 
-  int Counter = Count();
   if (ListOfLocalAtoms == NULL) { // allocated initial pointer
     // allocate and set each field to NULL
-    ListOfLocalAtoms = Malloc<atom**>(Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
-    if (ListOfLocalAtoms != NULL) {
-      for (int i = Counter; i--;)
-        ListOfLocalAtoms[i] = NULL;
-      FreeList = FreeList && true;
-    } else
+    const int Counter = Count();
+    ListOfLocalAtoms = Calloc<atom**>(Counter, "MoleculeLeafClass::FillListOfLocalAtoms - ***ListOfLocalAtoms");
+    if (ListOfLocalAtoms == NULL) {
+      FreeList = FreeList && false;
       status = false;
+    }
   }
 
@@ -961 +1040 @@
   if (FragmentList == NULL) {
     KeySetCounter = Count();
-    FragmentList = Malloc<Graph*>(KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
-    for (int i = KeySetCounter; i--;)
-      FragmentList[i] = NULL;
+    FragmentList = Calloc<Graph*>(KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
     KeySetCounter = 0;
   }