Changeset 183f35

Timestamp:

Apr 28, 2008, 9:22:38 AM (17 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:

b0a0c3

Parents:

08f1a0

Message:

FragmentMolecule(): Adjacency store/check and parsing of KeySetFile incorporated

• StoreFragmentFromKey_S_et() renamed
• StoreFragmentFromKeySet() now does not get whole config structure as parameter but only needed IsAngstroem value
• two new functions: ScanBufferIntoKeySet() - scans a text line of numbers into a KeySet deque ParseKeySetFile() - parses a file of the above lines, fills each into a KeySet and spills out fragment molecules by calling StoreFragmentFromKeySet()
• FragmentMolecule(): Adjacency file written, read and checked if it matches the just created bond structure. If so, then we don't fragment again but parse the KeySetFile and reconstruct fragments thereby.

Location:

src

Files:

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

src/molecules.cpp

@@ -1851 +1851 @@
 };  
 
+/** Scans a single line for number and puts them into \a KeySet.
+ * \param *out output stream for debugging
+ * \param *buffer buffer to scan
+ * \param &CurrentSet filled KeySet on return 
+ * \return true - at least one valid atom id parsed, false - CurrentSet is empty
+ */
+bool molecule::ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
+{
+  stringstream line;
+  int AtomNr;
+  int status = 0;
+  
+  line.str(buffer);
+  while (!line.eof()) {
+    line >> AtomNr;
+    if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+      CurrentSet.insert(AtomNr);
+      status++;
+    } // else it's "-1" or else and thus must not be added
+  }
+  return (status != 0);
+};
+
+/** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
+ * \param *out output stream for debugging
+ * \param *filename name of KeySet file
+ * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
+ * \param *FragmentList NULL, filled on return
+ * \param IsAngstroem whether we have Ansgtroem or bohrradius
+ * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
+ */
+bool molecule::ParseKeySetFile(ofstream *out, char *filename, atom **ListOfAtoms, MoleculeListClass *&FragmentList, bool IsAngstroem)
+{
+  bool status = true;
+  ifstream KeySetFile;
+  stringstream line;
+  
+  if (FragmentList != NULL) { // check list pointer
+    cerr << "Error: FragmentList was not NULL as supposed to be, already atoms present therein?" << endl;
+    return false;
+  }
+  cout << Verbose(1) << "Parsing the KeySet file ... ";
+  // open file and read
+  KeySetFile.open(filename);
+  if (KeySetFile != NULL) {
+    // each line represents a new fragment
+    int NumberOfFragments = 0;
+    char *buffer = (char *) Malloc(sizeof(char)*255, "molecule::ParseKeySetFile - *buffer");
+    // 1. scan through file to know number of fragments
+    while (!KeySetFile.eof()) {
+      KeySetFile.getline(buffer, 255);
+      if (strlen(buffer) > 0) // there is at least on possible number on the parsed line
+        NumberOfFragments++;
+    }
+    // 2. allocate the MoleculeListClass accordingly
+    FragmentList = new MoleculeListClass(NumberOfFragments, AtomCount);
+    // 3. Clear File, go to beginning and parse again, now adding found ids to each keyset and converting into molecules
+    KeySetFile.clear();
+    KeySetFile.seekg(ios::beg);
+    NumberOfFragments = 0;
+    while ((!KeySetFile.eof()) && (FragmentList->NumberOfMolecules > NumberOfFragments)) {
+      KeySetFile.getline(buffer, 255);
+      KeySet CurrentSet;
+      if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet)))  // if at least one valid atom was added, write config
+        FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
+    }
+    // 4. Free and done
+    Free((void **)&buffer, "molecule::ParseKeySetFile - *buffer");
+    cout << "done." << endl;
+  } else {
+    cout << "File not found." << endl;
+    status = false;
+  }
+  
+  return status;
+};
+
 /** Performs a many-body bond order analysis for a given bond order.
  * Writes for each fragment a config file.
@@ -1873 +1950 @@
   MoleculeLeafClass *MolecularWalker = NULL;
   MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
-
+  fstream File;
+  char *filename = NULL;
+  bool status = true;
+
+  cout << endl;
 #ifdef ADDHYDROGEN
   cout << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
@@ -1880 +1961 @@
 #endif
 
+  // fill the adjacency list
   CreateListOfBondsPerAtom(out);
 
-  // first perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs
-  Subgraphs = DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
-  MolecularWalker = Subgraphs;
-  // fill the bond structure of the individually stored subgraphs
-  while (MolecularWalker->next != NULL) {
-    MolecularWalker = MolecularWalker->next;
-    cout << Verbose(1) << "Creating adjacency list for subgraph " << MolecularWalker << "." << endl;
-    MolecularWalker->Leaf->CreateAdjacencyList((ofstream *)&cout, BondDistance);
-    MolecularWalker->Leaf->CreateListOfBondsPerAtom((ofstream *)&cout);
-  }
-  // fragment all subgraphs
-  if ((MinimumRingSize != -1) && ((BottomUpOrder >= MinimumRingSize) || (TopDownOrder >= MinimumRingSize))) {
-    cout << Verbose(0) << "Bond order greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
+  // === open and compare adjacency file ===
+  filename = (char *) Malloc(sizeof(char)*255, "molecule::FragmentMolecule - filename");
+  sprintf(filename, "%s/%s%s", configuration->GetDefaultPath(), "BondFragment", "Adjacency.dat");
+  File.open(filename, ios::in);
+  cout << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ...";
+  if (File != NULL) {
+    // allocate storage structure
+    atom **ListOfAtoms = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::FragmentMolecule - ListOfAtoms"); // will contain the atom ids whose bond partners differ
+    int NonMatchNumber = 0;   // will number of atoms with differing bond structure
+    int *CurrentBonds = (int *) Malloc(sizeof(int)*8, "molecule::FragmentMolecule - CurrentBonds"); // contains parsed bonds of current atom
+    int CurrentBondsOfAtom;
+    Walker = start;
+    while (Walker->next != end) {
+      Walker = Walker->next;
+      if ((Walker->nr >= 0) && (Walker->nr < AtomCount)) {
+        ListOfAtoms[Walker->nr] = Walker;
+      } else
+        break;
+    }
+    if (Walker->next == end) {  // everything went alright
+      // Parse the file line by line and count the bonds
+      while (!File.eof()) {
+        File.getline(filename, 255);
+        stringstream line;
+        line.str(filename);
+        int AtomNr = -1;
+        line >> AtomNr;
+        CurrentBondsOfAtom = -1; // we count one too far due to line end
+        // parse into structure
+        if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+          while (!line.eof())
+            line >> CurrentBonds[ ++CurrentBondsOfAtom ];
+          // compare against present bonds
+          Walker = ListOfAtoms[AtomNr];
+          //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
+          if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
+            for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
+              int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(Walker)->nr;
+              int j = 0;
+              for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
+              if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
+                ListOfAtoms[AtomNr] = NULL;
+                NonMatchNumber++;
+                status = false;
+                //cout << "[" << id << "]\t";
+              } else {
+                //cout << id << "\t";
+              }
+            }
+            //cout << endl;
+          } else {
+            cout << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
+            status = false;
+          }
+        }
+      }
+      File.close();
+      File.clear();
+      if (status) { // if equal we parse the KeySetFile
+        cout << " done: Equal." << endl;
+        sprintf(filename, "%s/%s%s", configuration->GetDefaultPath(), "BondFragment", "KeySets.dat");
+        status = ParseKeySetFile(out, filename, ListOfAtoms, FragmentList, configuration->GetIsAngstroem());
+      } else
+        cout << " done: Not equal by " << NonMatchNumber << " atoms." << endl;
+    } else {
+      cout << " range of nuclear ids exceeded [0, AtomCount)." << endl;
+      status = false;
+    }
+    Free((void **)&CurrentBonds, "molecule::FragmentMolecule - **CurrentBonds");
+    Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - *ListOfAtoms");
   } else {
-    FragmentCounter = 0;
+    cout << " Adjacency file not found." << endl;
+    status = false;
+  }
+
+  // =================================== Begin of FRAGMENTATION ===============================
+  if (!status) {
+    // === store Adjacency file ===
+    sprintf(filename, "%s/%s%s", configuration->GetDefaultPath(), "BondFragment", "Adjacency.dat");
+    File.open(filename, ios::out);
+    cout << Verbose(1) << "Saving adjacency list ... ";
+    if (File != NULL) {
+      Walker = start;
+      while(Walker->next != end) {
+        Walker = Walker->next;
+        File << Walker->nr << "\t";
+        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
+          File << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
+        File << endl;
+      }
+      File.close();
+      cout << "done." << endl;
+    } else {
+      cout << "failed." << endl;
+    }
+    Free((void **)&filename, "molecule::FragmentMolecule - filename");
+  
+    // === first perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs ===
+    Subgraphs = DepthFirstSearchAnalysis((ofstream *)&cout, false, MinimumRingSize);
     MolecularWalker = Subgraphs;
-    // count subgraphs
+    // fill the bond structure of the individually stored subgraphs
     while (MolecularWalker->next != NULL) {
       MolecularWalker = MolecularWalker->next;
-      FragmentCounter++;
-    }
-    BondFragments = (MoleculeListClass **) Malloc(sizeof(MoleculeListClass *)*FragmentCounter, "molecule::FragmentMolecule - **BondFragments");
-    // fill the bond fragment list
-    FragmentCounter = 0;
-    MolecularWalker = Subgraphs;
-    while (MolecularWalker->next != NULL) {
-      MolecularWalker = MolecularWalker->next;
-      cout << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
-      if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
-          // output ListOfBondsPerAtom for debugging
-        *out << Verbose(0) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
-        Walker = MolecularWalker->Leaf->start;
-        while (Walker->next != MolecularWalker->Leaf->end) {
-          Walker = Walker->next; 
-#ifdef ADDHYDROGEN
-          if (Walker->type->Z != 1) {   // regard only non-hydrogen
-#endif
-            *out << Verbose(0) << "Atom " << Walker->Name << " has Bonds: "<<endl;
-            for(int j=0;j<MolecularWalker->Leaf->NumberOfBondsPerAtom[Walker->nr];j++) {
-              *out << Verbose(1) << *(MolecularWalker->Leaf->ListOfBondsPerAtom)[Walker->nr][j] << endl;
+      cout << Verbose(1) << "Creating adjacency list for subgraph " << MolecularWalker << "." << endl;
+      MolecularWalker->Leaf->CreateAdjacencyList((ofstream *)&cout, BondDistance);
+      MolecularWalker->Leaf->CreateListOfBondsPerAtom((ofstream *)&cout);
+    }
+    
+    // === fragment all subgraphs ===
+    if ((MinimumRingSize != -1) && ((BottomUpOrder >= MinimumRingSize) || (TopDownOrder >= MinimumRingSize))) {
+      cout << Verbose(0) << "Bond order greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
+    } else {
+      FragmentCounter = 0;
+      MolecularWalker = Subgraphs;
+      // count subgraphs
+      while (MolecularWalker->next != NULL) {
+        MolecularWalker = MolecularWalker->next;
+        FragmentCounter++;
+      }
+      BondFragments = (MoleculeListClass **) Malloc(sizeof(MoleculeListClass *)*FragmentCounter, "molecule::FragmentMolecule - **BondFragments");
+      // fill the bond fragment list
+      FragmentCounter = 0;
+      MolecularWalker = Subgraphs;
+      while (MolecularWalker->next != NULL) {
+        MolecularWalker = MolecularWalker->next;
+        cout << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
+        if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
+            // output ListOfBondsPerAtom for debugging
+          *out << Verbose(0) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
+          Walker = MolecularWalker->Leaf->start;
+          while (Walker->next != MolecularWalker->Leaf->end) {
+            Walker = Walker->next; 
+  #ifdef ADDHYDROGEN
+            if (Walker->type->Z != 1) {   // regard only non-hydrogen
+  #endif
+              *out << Verbose(0) << "Atom " << Walker->Name << " has Bonds: "<<endl;
+              for(int j=0;j<MolecularWalker->Leaf->NumberOfBondsPerAtom[Walker->nr];j++) {
+                *out << Verbose(1) << *(MolecularWalker->Leaf->ListOfBondsPerAtom)[Walker->nr][j] << endl;
+              }
+  #ifdef ADDHYDROGEN
             }
-#ifdef ADDHYDROGEN
+  #endif
           }
-#endif
+          *out << endl;
+        
+          *out << Verbose(0) << endl << " ========== BOND ENERGY ========================= " << endl;
+          *out << Verbose(0) << "Begin of bond fragmentation." << endl;
+          BondFragments[FragmentCounter] = NULL;
+          if (Scheme == ANOVA) {
+            BondFragments[FragmentCounter] = MolecularWalker->Leaf->FragmentBOSSANOVA(out,BottomUpOrder,configuration);
+          }  
+          if ((Scheme == BottomUp) || (Scheme == Combined)) { // get overlapping subgraphs
+            BondFragments[FragmentCounter] = FragmentList = MolecularWalker->Leaf->FragmentBottomUp(out, BottomUpOrder, configuration, CutCyclic);
+          }
+          if (Scheme == TopDown) { // initialise top level with whole molecule
+            *out << Verbose(2) << "Initial memory allocating and initialising for whole molecule." << endl;
+            FragmentList = new MoleculeListClass(1, MolecularWalker->Leaf->AtomCount);
+            FragmentList->ListOfMolecules[0] = MolecularWalker->Leaf->CopyMolecule();
+            FragmentList->TEList[0] = 1.;
+          }
+          if ((Scheme == Combined) || (Scheme == TopDown)) {
+            *out << Verbose(1) << "Calling TopDown." << endl;
+            BondFragments[FragmentCounter] = FragmentList->FragmentTopDown(out, TopDownOrder, BondDistance, configuration, CutCyclic);
+            // remove this molecule from list again and free wrapper list
+            delete(FragmentList);
+            FragmentList = NULL;
+          }
+        } else {
+          cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
         }
-        *out << endl;
-      
-        *out << Verbose(0) << endl << " ========== BOND ENERGY ========================= " << endl;
-        *out << Verbose(0) << "Begin of bond fragmentation." << endl;
-        BondFragments[FragmentCounter] = NULL;
-        if (Scheme == ANOVA) {
-          BondFragments[FragmentCounter] = MolecularWalker->Leaf->FragmentBOSSANOVA(out,BottomUpOrder,configuration);
-        }  
-        if ((Scheme == BottomUp) || (Scheme == Combined)) { // get overlapping subgraphs
-          BondFragments[FragmentCounter] = FragmentList = MolecularWalker->Leaf->FragmentBottomUp(out, BottomUpOrder, configuration, CutCyclic);
-        }
-        if (Scheme == TopDown) { // initialise top level with whole molecule
-          *out << Verbose(2) << "Initial memory allocating and initialising for whole molecule." << endl;
-          FragmentList = new MoleculeListClass(1, MolecularWalker->Leaf->AtomCount);
-          FragmentList->ListOfMolecules[0] = MolecularWalker->Leaf->CopyMolecule();
-          FragmentList->TEList[0] = 1.;
-        }
-        if ((Scheme == Combined) || (Scheme == TopDown)) {
-          *out << Verbose(1) << "Calling TopDown." << endl;
-          BondFragments[FragmentCounter] = FragmentList->FragmentTopDown(out, TopDownOrder, BondDistance, configuration, CutCyclic);
-          // remove this molecule from list again and free wrapper list
-          delete(FragmentList);
-          FragmentList = NULL;
-        }
-      } else {
-        cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
-      }
-      TotalFragmentCounter += BondFragments[FragmentCounter]->NumberOfMolecules;
-      FragmentCounter++;  // next fragment list
-    }
-  }
-
-  // combine bond fragments list into a single one
-  FragmentList = new MoleculeListClass(TotalFragmentCounter, AtomCount);
-  TotalFragmentCounter = 0;
-  for (int i=0;i<FragmentCounter;i++) {
-    for(int j=0;j<BondFragments[i]->NumberOfMolecules;j++) {
-      FragmentList->ListOfMolecules[TotalFragmentCounter] =  BondFragments[i]->ListOfMolecules[j];
-      BondFragments[i]->ListOfMolecules[j] = NULL;
-      FragmentList->TEList[TotalFragmentCounter++] = BondFragments[i]->TEList[j];
-    }
-    delete(BondFragments[i]);
-  }
-  Free((void **)&BondFragments, "molecule::FragmentMolecule - **BondFragments");
-  
-  // now if there are actually any fragments to save on disk ...
+        TotalFragmentCounter += BondFragments[FragmentCounter]->NumberOfMolecules;
+        FragmentCounter++;  // next fragment list
+      }
+    }
+  
+    // === combine bond fragments list into a single one ===
+    FragmentList = new MoleculeListClass(TotalFragmentCounter, AtomCount);
+    TotalFragmentCounter = 0;
+    for (int i=0;i<FragmentCounter;i++) {
+      for(int j=0;j<BondFragments[i]->NumberOfMolecules;j++) {
+        FragmentList->ListOfMolecules[TotalFragmentCounter] =  BondFragments[i]->ListOfMolecules[j];
+        BondFragments[i]->ListOfMolecules[j] = NULL;
+        FragmentList->TEList[TotalFragmentCounter++] = BondFragments[i]->TEList[j];
+      }
+      delete(BondFragments[i]);
+    }
+    Free((void **)&BondFragments, "molecule::FragmentMolecule - **BondFragments");
+  } else
+    cout << Verbose(0) << "Using fragments reconstructed from the KeySetFile." << endl;
+  // ==================================== End of FRAGMENTATION ================================
+  
+  // === Save fragments' configuration to disk ===
   if (FragmentList != NULL) {
     // create a SortIndex to map from BFS labels to the sequence in which the atoms are given in the config file
@@ -2008 +2191 @@
     *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
   // free subgraph memory again
-  while (Subgraphs->next != NULL) {
-    Subgraphs = Subgraphs->next;
-    delete(Subgraphs->previous);
-  }
-  delete(Subgraphs);
+  if (Subgraphs != NULL) {
+    while (Subgraphs->next != NULL) {
+      Subgraphs = Subgraphs->next;
+      delete(Subgraphs->previous);
+    }
+    delete(Subgraphs);
+  }
 
   *out << Verbose(0) << "End of bond fragmentation." << endl;
@@ -2553 +2738 @@
 };
 
-/** Stores a fragment from \a SnakeStack into \a molecule.
+/** Stores a fragment from \a KeySet into \a molecule.
  * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
  * molecule and adds missing hydrogen where bonds were cut. 
  * \param *out output stream for debugging messages
  * \param &Leaflet pointer to KeySet structure 
- * \param *configuration configuration for writing config files for each fragment
+ * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return pointer to constructed molecule
  */
-molecule * molecule::StoreFragmentFromKeyset(ofstream *&out, KeySet &Leaflet, config *&configuration)
+molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
 {
   atom *Runner = NULL, *FatherOfRunner = NULL, *OtherFather = NULL;
@@ -2606 +2791 @@
 #ifdef ADDHYDROGEN
           *out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], configuration->GetIsAngstroem());
+          Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem);
 #endif
           //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
@@ -3557 +3742 @@
   for(Graph::iterator runner = FragmentList->begin(); runner != FragmentList->end(); runner++) {
     KeySet test = (*runner).first;
-    FragmentMoleculeList->ListOfMolecules[k] = StoreFragmentFromKeyset(out, test, configuration);
+    FragmentMoleculeList->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
     FragmentMoleculeList->TEList[k] = ((*runner).second).second;
     k++;

src/molecules.hpp

@@ -520 +520 @@
   /// Fragment molecule by two different approaches:
   void FragmentMolecule(ofstream *out, int BottomUpOrder, int TopDownOrder, enum BondOrderScheme Scheme, config *configuration, enum CutCyclicBond CutCyclic);
+  bool molecule::ParseKeySetFile(ofstream *out, char *filename, atom **ListOfAtoms, MoleculeListClass *&FragmentList, bool IsAngstroem);
+  bool molecule::ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
   MoleculeListClass * GetAtomicFragments(ofstream *out, int NumberOfTopAtoms, bool IsAngstroem, double factor, enum CutCyclicBond CutCyclic);
   void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem, enum CutCyclicBond CutCyclic);
@@ -531 +533 @@
         int CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, Graph *ListOfGraph, KeySet Fragment, double TEFactor, config *configuration);
   bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
-  molecule * StoreFragmentFromKeyset(ofstream *&out, KeySet &Leaflet, config *&configuration);
+  molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
   void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
   int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);