Changes in src/molecule_fragmentation.cpp [d4c9ae:35b698]

File:

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

src/molecule_fragmentation.cpp

@@ -5 +5 @@
  *      Author: heber
  */
+
+#include "Helpers/MemDebug.hpp"
 
 #include <cstring>
@@ -39 +41 @@
   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));
@@ -82 +82 @@
  * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
  * Finally, the temporary graph is inserted into the given \a FragmentList for return.
- * \param *out output stream for debugging
- * \param *path path to file
+ * \param &path path to file
  * \param *FragmentList empty, filled on return
  * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
  */
-bool ParseKeySetFile(char *path, Graph *&FragmentList)
+bool ParseKeySetFile(std::string &path, Graph *&FragmentList)
 {
   bool status = true;
@@ -94 +93 @@
   GraphTestPair testGraphInsert;
   int NumberOfFragments = 0;
-  char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
+  string filename;
 
   if (FragmentList == NULL) { // check list pointer
@@ -102 +101 @@
   // 1st pass: open file and read
   DoLog(1) && (Log() << Verbose(1) << "Parsing the KeySet file ... " << endl);
-  sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
-  InputFile.open(filename);
-  if (InputFile != NULL) {
+  filename = path + KEYSETFILE;
+  InputFile.open(filename.c_str());
+  if (InputFile.good()) {
     // each line represents a new fragment
-    char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - *buffer");
+    char buffer[MAXSTRINGSIZE];
     // 1. parse keysets and insert into temp. graph
     while (!InputFile.eof()) {
@@ -122 +121 @@
     InputFile.close();
     InputFile.clear();
-    Free(&buffer);
-    DoLog(1) && (Log() << Verbose(1) << "done." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "\t ... done." << endl);
   } else {
-    DoLog(1) && (Log() << Verbose(1) << "File " << filename << " not found." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "\t ... File " << filename << " not found." << endl);
     status = false;
   }
 
-  Free(&filename);
   return status;
 };
@@ -148 +145 @@
   int NumberOfFragments = 0;
   double TEFactor;
-  char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseTEFactorsFile - filename");
+  char filename[MAXSTRINGSIZE];
 
   if (FragmentList == NULL) { // check list pointer
@@ -179 +176 @@
   }
 
-  // free memory
-  Free(&filename);
-
   return status;
 };
 
 /** Stores key sets to file.
- * \param *out output stream for debugging
  * \param KeySetList Graph with Keysets
- * \param *path path to file
+ * \param &path path to file
  * \return true - file written successfully, false - writing failed
  */
-bool StoreKeySetFile(Graph &KeySetList, char *path)
-{
-  ofstream output;
+bool StoreKeySetFile(Graph &KeySetList, std::string &path)
+{
   bool status =  true;
-  string line;
+  string line = path + KEYSETFILE;
+  ofstream output(line.c_str());
 
   // open KeySet file
-  line = path;
-  line.append("/");
-  line += FRAGMENTPREFIX;
-  line += KEYSETFILE;
-  output.open(line.c_str(), ios::out);
   DoLog(1) && (Log() << Verbose(1) << "Saving key sets of the total graph ... ");
-  if(output != NULL) {
+  if(output.good()) {
     for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
       for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
@@ -307 +295 @@
 
 /** Scans the adaptive order file and insert (index, value) into map.
- * \param *out output stream for debugging
- * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
+ * \param &path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
  * \param &IndexedKeySetList list to find key set for a given index \a No
  * \return adaptive criteria list from file
  */
-map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(char *path, map<int,KeySet> &IndexKeySetList)
+map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(std::string &path, map<int,KeySet> &IndexKeySetList)
 {
   map<int, pair<double,int> > *AdaptiveCriteriaList = new map<int, pair<double,int> >;
   int No = 0, FragOrder = 0;
   double Value = 0.;
-  char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
-  sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
-  ifstream InputFile(buffer, ios::in);
+  char buffer[MAXSTRINGSIZE];
+  string filename = path + ENERGYPERFRAGMENT;
+  ifstream InputFile(filename.c_str());
+
+  if (InputFile.fail()) {
+    DoeLog(1) && (eLog() << Verbose(1) << "Cannot find file " << filename << "." << endl);
+    return AdaptiveCriteriaList;
+  }
 
   if (CountLinesinFile(InputFile) > 0) {
@@ -345 +337 @@
     InputFile.clear();
   }
-  Free(&buffer);
 
   return AdaptiveCriteriaList;
@@ -359 +350 @@
 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);
@@ -391 +382 @@
 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;
@@ -425 +416 @@
 
 /** Checks whether the OrderAtSite is still below \a Order at some site.
- * \param *out output stream for debugging
  * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
  * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
  * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
  * \param *MinimumRingSize array of max. possible order to avoid loops
- * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
+ * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
  * \return true - needs further fragmentation, false - does not need fragmentation
  */
-bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
-{
-  atom *Walker = start;
+bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, std::string path)
+{
   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;
 
@@ -457 +446 @@
     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;
         }
@@ -474 +462 @@
     MarkUpdateCandidates(AtomMask, *FinalRootCandidates, Order, this);
 
-    Free(&IndexKeySetList);
-    Free(&AdaptiveCriteriaList);
-    Free(&FinalRootCandidates);
+    delete[](IndexKeySetList);
+    delete[](AdaptiveCriteriaList);
+    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;
 
@@ -500 +487 @@
   }
 
-  PrintAtomMask(AtomMask, AtomCount); // for debugging
+  PrintAtomMask(AtomMask, getAtomCount()); // for debugging
 
   return status;
@@ -516 +503 @@
     return false;
   }
-  SortIndex = Malloc<int>(AtomCount, "molecule::CreateMappingLabelsToConfigSequence: *SortIndex");
-  for(int i=AtomCount;i--;)
+  SortIndex = new int[getAtomCount()];
+  for(int i=getAtomCount();i--;)
     SortIndex[i] = -1;
 
@@ -524 +511 @@
 
   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;
 };
 
@@ -539 +581 @@
  * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
  * subgraph in the MoleculeListClass.
- * \param *out output stream for debugging
  * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
- * \param *configuration configuration for writing config files for each fragment
+ * \param &prefix path and prefix of the bond order configs to be written
  * \return 1 - continue, 2 - stop (no fragmentation occured)
  */
-int molecule::FragmentMolecule(int Order, config *configuration)
+int molecule::FragmentMolecule(int Order, std::string &prefix)
 {
   MoleculeListClass *BondFragments = NULL;
-  int *SortIndex = NULL;
-  int *MinimumRingSize = new int[AtomCount];
+  int *MinimumRingSize = new int[getAtomCount()];
   int FragmentCounter;
   MoleculeLeafClass *MolecularWalker = NULL;
@@ -576 +616 @@
 
   // create lookup table for Atom::nr
-  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(start, end, ListOfAtoms, AtomCount);
+  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, getAtomCount());
 
   // === compare it with adjacency file ===
-  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(configuration->configpath, ListOfAtoms);
-  Free(&ListOfAtoms);
+  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(prefix, ListOfAtoms);
+  delete[](ListOfAtoms);
 
   // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
@@ -586 +626 @@
 
   // 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;
@@ -598 +638 @@
 //    // 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";
@@ -613 +653 @@
 
   // ===== 3. if structure still valid, parse key set file and others =====
-  FragmentationToDo = FragmentationToDo && ParseKeySetFile(configuration->configpath, ParsedFragmentList);
+  FragmentationToDo = FragmentationToDo && ParseKeySetFile(prefix, ParsedFragmentList);
 
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
-  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(configuration->configpath);
+  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(prefix);
 
   // =================================== Begin of FRAGMENTATION ===============================
@@ -624 +664 @@
   // ===== 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))) {
+  while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, prefix))) {
     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));
@@ -640 +680 @@
       DoLog(1) && (Log() << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl);
       //MolecularWalker->Leaf->OutputListOfBonds(out);  // output atom::ListOfBonds for debugging
-      if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
+      if (MolecularWalker->Leaf->hasBondStructure()) {
         // call BOSSANOVA method
         DoLog(0) && (Log() << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl);
@@ -672 +712 @@
     delete(Subgraphs);
   }
-  Free(&FragmentList);
+  delete[](FragmentList);
 
   // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
@@ -682 +722 @@
     KeySet test = (*runner).first;
     DoLog(0) && (Log() << Verbose(0) << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl);
-    BondFragments->insert(StoreFragmentFromKeySet(test, configuration));
+    BondFragments->insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
     k++;
   }
@@ -690 +730 @@
   if (BondFragments->ListOfMolecules.size() != 0) {
     // create the SortIndex from BFS labels to order in the config file
+    int *SortIndex = NULL;
     CreateMappingLabelsToConfigSequence(SortIndex);
 
     DoLog(1) && (Log() << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl);
-    if (BondFragments->OutputConfigForListOfFragments(configuration, SortIndex))
+    if (BondFragments->OutputConfigForListOfFragments(prefix, SortIndex))
       DoLog(1) && (Log() << Verbose(1) << "All configs written." << endl);
     else
@@ -699 +740 @@
 
     // store force index reference file
-    BondFragments->StoreForcesFile(configuration->configpath, SortIndex);
+    BondFragments->StoreForcesFile(prefix, SortIndex);
 
     // store keysets file
-    StoreKeySetFile(TotalGraph, configuration->configpath);
-
-    // store Adjacency file
-    char *filename = Malloc<char> (MAXSTRINGSIZE, "molecule::FragmentMolecule - *filename");
-    strcpy(filename, FRAGMENTPREFIX);
-    strcat(filename, ADJACENCYFILE);
-    StoreAdjacencyToFile(configuration->configpath, filename);
-    Free(&filename);
+    StoreKeySetFile(TotalGraph, prefix);
+
+    {
+      // store Adjacency file
+      std::string filename = prefix + ADJACENCYFILE;
+      StoreAdjacencyToFile(filename);
+    }
 
     // store Hydrogen saturation correction file
-    BondFragments->AddHydrogenCorrection(configuration->configpath);
+    BondFragments->AddHydrogenCorrection(prefix);
 
     // store adaptive orders into file
-    StoreOrderAtSiteFile(configuration->configpath);
+    StoreOrderAtSiteFile(prefix);
 
     // restore orbital and Stop values
-    CalculateOrbitals(*configuration);
+    //CalculateOrbitals(*configuration);
 
     // free memory for bond part
     DoLog(1) && (Log() << Verbose(1) << "Freeing bond memory" << endl);
-    Free(&FragmentList); // remove bond molecule from memory
-    Free(&SortIndex);
+    delete[](SortIndex);
   } else {
     DoLog(1) && (Log() << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl);
@@ -736 +775 @@
 /** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
  * Atoms not present in the file get "-1".
- * \param *out output stream for debugging
- * \param *path path to file ORDERATSITEFILE
+ * \param &path path to file ORDERATSITEFILE
  * \return true - file writable, false - not writable
  */
-bool molecule::StoreOrderAtSiteFile(char *path)
-{
-  stringstream line;
+bool molecule::StoreOrderAtSiteFile(std::string &path)
+{
+  string line;
   ofstream file;
 
-  line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
-  file.open(line.str().c_str());
+  line = path + ORDERATSITEFILE;
+  file.open(line.c_str());
   DoLog(1) && (Log() << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl);
-  if (file != NULL) {
+  if (file.good()) {
     ActOnAllAtoms( &atom::OutputOrder, &file );
     file.close();
@@ -754 +792 @@
     return true;
   } else {
-    DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line.str() << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line << "." << endl);
     return false;
   }
@@ -761 +799 @@
 /** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
  * Atoms not present in the file get "0".
- * \param *out output stream for debugging
- * \param *path path to file ORDERATSITEFILEe
+ * \param &path path to file ORDERATSITEFILEe
  * \return true - file found and scanned, false - file not found
  * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
  */
-bool molecule::ParseOrderAtSiteFromFile(char *path)
-{
-  unsigned char *OrderArray = Calloc<unsigned char>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
-  bool *MaxArray = Calloc<bool>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
+bool molecule::ParseOrderAtSiteFromFile(std::string &path)
+{
+  unsigned char *OrderArray = new unsigned char[getAtomCount()];
+  bool *MaxArray = new bool[getAtomCount()];
   bool status;
   int AtomNr, value;
-  stringstream line;
+  string line;
   ifstream file;
 
+  for(int i=0;i<getAtomCount();i++) {
+    OrderArray[i] = 0;
+    MaxArray[i] = false;
+  }
+
   DoLog(1) && (Log() << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl);
-  line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
-  file.open(line.str().c_str());
-  if (file != NULL) {
+  line = path + ORDERATSITEFILE;
+  file.open(line.c_str());
+  if (file.good()) {
     while (!file.eof()) { // parse from file
       AtomNr = -1;
@@ -796 +838 @@
     SetAtomValueToIndexedArray( MaxArray, &atom::nr, &atom::MaxOrder );
 
-    DoLog(1) && (Log() << Verbose(1) << "done." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "\t ... done." << endl);
     status = true;
   } else {
-    DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line.str() << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "\t ... failed to open file " << line << "." << endl);
     status = false;
   }
-  Free(&OrderArray);
-  Free(&MaxArray);
+  delete[](OrderArray);
+  delete[](MaxArray);
 
   DoLog(1) && (Log() << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl);
@@ -872 +914 @@
   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
@@ -889 +933 @@
 //            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;
@@ -899 +943 @@
 //          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
   }
@@ -929 +975 @@
 molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
 {
-  atom **SonList = Calloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList");
+  atom **SonList = new atom*[getAtomCount()];
   molecule *Leaf = World::getInstance().createMolecule();
+
+  for(int i=0;i<getAtomCount();i++)
+    SonList[i] = NULL;
 
 //  Log() << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
@@ -939 +988 @@
 
   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
-  Free(&SonList);
+  delete[](SonList);
 //  Log() << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
   return Leaf;
@@ -1084 +1133 @@
   int bits, TouchedIndex, SubSetDimension, SP, Added;
   int SpaceLeft;
-  int *TouchedList = Malloc<int>(SubOrder + 1, "molecule::SPFragmentGenerator: *TouchedList");
-  bond **BondsList = NULL;
+  int *TouchedList = new int[SubOrder + 1];
   KeySetTestPair TestKeySetInsert;
 
@@ -1124 +1172 @@
 
           // then allocate and fill the list
-          BondsList = Malloc<bond*>(SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
+          bond *BondsList[SubSetDimension];
           SubSetDimension = FillBondsList(BondsList, FragmentSearch->BondsPerSPList[2*SP], FragmentSearch->BondsPerSPList[2*SP+1], TouchedList, TouchedIndex);
 
@@ -1130 +1178 @@
           Log() << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
           SPFragmentGenerator(FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
-
-          Free(&BondsList);
         }
       } else {
@@ -1153 +1199 @@
     }
   }
-  Free(&TouchedList);
+  delete[](TouchedList);
   Log() << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
 };
@@ -1165 +1211 @@
 void InitialiseSPList(int Order, struct UniqueFragments &FragmentSearch)
 {
-  FragmentSearch.BondsPerSPList = Malloc<bond*>(Order * 2, "molecule::PowerSetGenerator: ***BondsPerSPList");
-  FragmentSearch.BondsPerSPCount = Malloc<int>(Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
+  FragmentSearch.BondsPerSPList = new bond* [Order * 2];
+  FragmentSearch.BondsPerSPCount = new int[Order];
   for (int i=Order;i--;) {
     FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
@@ -1184 +1230 @@
 void FreeSPList(int Order, struct UniqueFragments &FragmentSearch)
 {
-  Free(&FragmentSearch.BondsPerSPCount);
+  delete[](FragmentSearch.BondsPerSPCount);
   for (int i=Order;i--;) {
     delete(FragmentSearch.BondsPerSPList[2*i]);
     delete(FragmentSearch.BondsPerSPList[2*i+1]);
   }
-  Free(&FragmentSearch.BondsPerSPList);
+  delete[](FragmentSearch.BondsPerSPList);
 };
 
@@ -1370 +1416 @@
 int molecule::PowerSetGenerator(int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
 {
-  bond **BondsList = NULL;
   int Counter = FragmentSearch.FragmentCounter; // mark current value of counter
 
@@ -1394 +1439 @@
 
     // prepare the subset and call the generator
-    BondsList = Calloc<bond*>(FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
+    bond* BondsList[FragmentSearch.BondsPerSPCount[0]];
+    for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++)
+      BondsList[i] = NULL;
     BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond
 
     SPFragmentGenerator(&FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
-
-    Free(&BondsList);
   } else {
     DoLog(0) && (Log() << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl);
@@ -1507 +1552 @@
       }
     }
-    Free(&FragmentLowerOrdersList[RootNr]);
+    delete[](FragmentLowerOrdersList[RootNr]);
     RootNr++;
   }
-  Free(&FragmentLowerOrdersList);
+  delete[](FragmentLowerOrdersList);
 };
 
@@ -1549 +1594 @@
   // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
   // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
-  NumMoleculesOfOrder = Calloc<int>(UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
-  FragmentLowerOrdersList = Calloc<Graph**>(UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
+  NumMoleculesOfOrder = new int[UpgradeCount];
+  FragmentLowerOrdersList = new Graph**[UpgradeCount];
+
+  for(int i=0;i<UpgradeCount;i++) {
+    NumMoleculesOfOrder[i] = 0;
+    FragmentLowerOrdersList[i] = NULL;
+  }
 
   // initialise the fragments structure
@@ -1556 +1606 @@
   FragmentSearch.FragmentSet = new KeySet;
   FragmentSearch.Root = FindAtom(RootKeyNr);
-  FragmentSearch.ShortestPathList = Malloc<int>(AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
-  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);
   }
 
@@ -1577 +1625 @@
     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])) {
@@ -1592 +1640 @@
       // allocate memory for all lower level orders in this 1D-array of ptrs
       NumLevels = 1 << (Order-1); // (int)pow(2,Order);
-      FragmentLowerOrdersList[RootNr] = Calloc<Graph*>(NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
+      FragmentLowerOrdersList[RootNr] = new Graph*[NumLevels];
+      for (int i=0;i<NumLevels;i++)
+        FragmentLowerOrdersList[RootNr][i] = NULL;
 
       // create top order where nothing is reduced
@@ -1628 +1678 @@
 
   // cleanup FragmentSearch structure
-  Free(&FragmentSearch.ShortestPathList);
+  delete[](FragmentSearch.ShortestPathList);
   delete(FragmentSearch.FragmentSet);
 
@@ -1641 +1691 @@
   CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, this);
   FreeAllOrdersList(FragmentLowerOrdersList, RootStack, this);
-  Free(&NumMoleculesOfOrder);
+  delete[](NumMoleculesOfOrder);
 
   DoLog(0) && (Log() << Verbose(0) << "End of FragmentBOSSANOVA." << endl);
@@ -1664 +1714 @@
   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 = Calloc<enum Shading>(AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");
+  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
     Translationvector.Zero();
     // scan all bonds
-    Binder = first;
     flag = false;
-    while ((!flag) && (Binder->next != last)) {
-      Binder = Binder->next;
-      for (int i=NDIM;i--;) {
-        tmp = fabs(Binder->leftatom->x[i] - Binder->rightatom->x[i]);
-        //Log() << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
-        if (tmp > BondDistance) {
-          OtherBinder = Binder->next; // note down binding partner for later re-insertion
-          unlink(Binder);   // unlink bond
-          DoLog(2) && (Log() << Verbose(2) << "Correcting at bond " << *Binder << "." << endl);
-          flag = true;
-          break;
+    for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner)
+      for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); (!flag) && (BondRunner != (*AtomRunner)->ListOfBonds.end()); ++BondRunner) {
+        Binder = (*BondRunner);
+        for (int i=NDIM;i--;) {
+          tmp = fabs(Binder->leftatom->x[i] - Binder->rightatom->x[i]);
+          //Log() << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
+          if (tmp > BondDistance) {
+            OtherBinder = Binder->next; // note down binding partner for later re-insertion
+            unlink(Binder);   // unlink bond
+            DoLog(2) && (Log() << Verbose(2) << "Correcting at bond " << *Binder << "." << endl);
+            flag = true;
+            break;
+          }
         }
       }
-    }
     if (flag) {
       // create translation vector from their periodically modified distance
@@ -1701 +1753 @@
       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);
@@ -1727 +1779 @@
   // free allocated space from ReturnFullMatrixforSymmetric()
   delete(AtomStack);
-  Free(&ColorList);
-  Free(&matrix);
+  delete[](ColorList);
+  delete[](matrix);
   DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
 };