Changes in src/moleculelist.cpp [2746be:437922]

File:

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

src/moleculelist.cpp

@@ -1 +1 @@
 /** \file MoleculeListClass.cpp
- * 
+ *
  * Function implementations for the class MoleculeListClass.
- * 
+ *
  */
 
@@ -13 +13 @@
 MoleculeListClass::MoleculeListClass()
 {
-};
-
-/** constructor for MoleculeListClass.
- * \param NumMolecules number of molecules to allocate for
- * \param NumAtoms number of atoms to allocate for
- */
-MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
-{
-  ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
-  for (int i=NumMolecules;i--;)
-    ListOfMolecules[i] = NULL;
-  NumberOfMolecules = NumMolecules;
-  NumberOfTopAtoms = NumAtoms;
-};
-
+  // empty lists
+  ListOfMolecules.clear();
+  MaxIndex = 1;
+};
 
 /** Destructor for MoleculeListClass.
@@ -34 +23 @@
 {
   cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
-  for (int i=NumberOfMolecules;i--;) {
-    if (ListOfMolecules[i] != NULL) { // if NULL don't free
-      cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
-      delete(ListOfMolecules[i]);
-      ListOfMolecules[i] = NULL;
-    }
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+    cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl;
+    delete (*ListRunner);
   }
   cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
-  Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
+  ListOfMolecules.clear(); // empty list
+};
+
+/** Insert a new molecule into the list and set its number.
+ * \param *mol molecule to add to list.
+ * \return true - add successful
+ */
+void MoleculeListClass::insert(molecule *mol)
+{
+  mol->IndexNr = MaxIndex++;
+  ListOfMolecules.push_back(mol);
 };
 
@@ -57 +53 @@
   atom *aWalker = NULL;
   atom *bWalker = NULL;
-  
+
   // sort each atom list and put the numbers into a list, then go through
   //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
-  if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
+  if ((**(molecule **) a).AtomCount < (**(molecule **) b).AtomCount) {
     return -1;
-  } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
-    return +1;
+  } else {
+    if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount)
+      return +1;
     else {
-      Count = (**(molecule **)a).AtomCount;
+      Count = (**(molecule **) a).AtomCount;
       aList = new int[Count];
       bList = new int[Count];
-  
+
       // fill the lists
-      aWalker = (**(molecule **)a).start;
-      bWalker = (**(molecule **)b).start;
+      aWalker = (**(molecule **) a).start;
+      bWalker = (**(molecule **) b).start;
       Counter = 0;
       aCounter = 0;
       bCounter = 0;
-      while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
+      while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
         aWalker = aWalker->next;
         bWalker = bWalker->next;
@@ -90 +87 @@
       // check if AtomCount was for real
       flag = 0;
-      if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
+      if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
         flag = -1;
       } else {
-        if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
+        if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end))
           flag = 1;
       }
       if (flag == 0) {
         // sort the lists
-        gsl_heapsort(aList,Count, sizeof(int), CompareDoubles); 
-        gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
-        // compare the lists 
-        
+        gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
+        gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
+        // compare the lists
+
         flag = 0;
-        for(int i=0;i<Count;i++) {
+        for (int i = 0; i < Count; i++) {
           if (aList[i] < bList[i]) {
             flag = -1;
@@ -114 +111 @@
         }
       }
-      delete[](aList);
-      delete[](bList);
+      delete[] (aList);
+      delete[] (bList);
       return flag;
     }
   }
-  return  -1;
+  return -1;
+};
+
+/** Output of a list of all molecules.
+ * \param *out output stream
+ */
+void MoleculeListClass::Enumerate(ofstream *out)
+{
+  int i=1;
+  element* Elemental = NULL;
+  atom *Walker = NULL;
+  int Counts[MAX_ELEMENTS];
+  double size=0;
+  Vector Origin;
+
+  // header
+  *out << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
+  cout << Verbose(0) << "-----------------------------------------------" << endl;
+  if (ListOfMolecules.size() == 0)
+    *out << "\tNone" << endl;
+  else {
+    Origin.Zero();
+    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+      // reset element counts
+      for (int j = 0; j<MAX_ELEMENTS;j++)
+        Counts[j] = 0;
+      // count atoms per element and determine size of bounding sphere
+      size=0.;
+      Walker = (*ListRunner)->start;
+      while (Walker->next != (*ListRunner)->end) {
+        Walker = Walker->next;
+        Counts[Walker->type->Z]++;
+        if (Walker->x.DistanceSquared(&Origin) > size)
+          size = Walker->x.DistanceSquared(&Origin);
+      }
+      // output Index, Name, number of atoms, chemical formula
+      *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
+      Elemental = (*ListRunner)->elemente->end;
+      while(Elemental->previous != (*ListRunner)->elemente->start) {
+        Elemental = Elemental->previous;
+        if (Counts[Elemental->Z] != 0)
+          *out << Elemental->symbol << Counts[Elemental->Z];
+      }
+      // Center and size
+      *out << "\t" << (*ListRunner)->Center << "\t" << sqrt(size) << endl;
+    }
+  }
+};
+
+/** Returns the molecule with the given index \a index.
+ * \param index index of the desired molecule
+ * \return pointer to molecule structure, NULL if not found
+ */
+molecule * MoleculeListClass::ReturnIndex(int index)
+{
+  for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
+    if ((*ListRunner)->IndexNr == index)
+      return (*ListRunner);
+  return NULL;
+};
+
+/** Simple merge of two molecules into one.
+ * \param *mol destination molecule
+ * \param *srcmol source molecule
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ */
+bool MoleculeListClass::SimpleMerge(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+    return false;
+
+  // put all molecules of src into mol
+  atom *Walker = srcmol->start;
+  atom *NextAtom = Walker->next;
+  while (NextAtom != srcmol->end) {
+    Walker = NextAtom;
+    NextAtom = Walker->next;
+    srcmol->UnlinkAtom(Walker);
+    mol->AddAtom(Walker);
+  }
+
+  // remove src
+  ListOfMolecules.remove(srcmol);
+  delete(srcmol);
+  return true;
+};
+
+/** Simple add of one molecules into another.
+ * \param *mol destination molecule
+ * \param *srcmol source molecule
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ */
+bool MoleculeListClass::SimpleAdd(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+    return false;
+
+  // put all molecules of src into mol
+  atom *Walker = srcmol->start;
+  atom *NextAtom = Walker->next;
+  while (NextAtom != srcmol->end) {
+    Walker = NextAtom;
+    NextAtom = Walker->next;
+    Walker = mol->AddCopyAtom(Walker);
+    Walker->father = Walker;
+  }
+
+  return true;
+};
+
+/** Simple merge of a given set of molecules into one.
+ * \param *mol destination molecule
+ * \param *src index of set of source molecule
+ * \param N number of source molecules
+ * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
+ */
+bool MoleculeListClass::SimpleMultiMerge(molecule *mol, int *src, int N)
+{
+  bool status = true;
+  // check presence of all source molecules
+  for (int i=0;i<N;i++) {
+    molecule *srcmol = ReturnIndex(src[i]);
+    status = status && SimpleMerge(mol, srcmol);
+  }
+  return status;
+};
+
+/** Simple add of a given set of molecules into one.
+ * \param *mol destination molecule
+ * \param *src index of set of source molecule
+ * \param N number of source molecules
+ * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
+ */
+bool MoleculeListClass::SimpleMultiAdd(molecule *mol, int *src, int N)
+{
+  bool status = true;
+  // check presence of all source molecules
+  for (int i=0;i<N;i++) {
+    molecule *srcmol = ReturnIndex(src[i]);
+    status = status && SimpleAdd(mol, srcmol);
+  }
+  return status;
+};
+
+/** Scatter merge of a given set of molecules into one.
+ * Scatter merge distributes the molecules in such a manner that they don't overlap.
+ * \param *mol destination molecule
+ * \param *src index of set of source molecule
+ * \param N number of source molecules
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ * \TODO find scatter center for each src molecule
+ */
+bool MoleculeListClass::ScatterMerge(molecule *mol, int *src, int N)
+{
+  // check presence of all source molecules
+  for (int i=0;i<N;i++) {
+    // get pointer to src molecule
+    molecule *srcmol = ReturnIndex(src[i]);
+    if (srcmol == NULL)
+      return false;
+  }
+  // adapt each Center
+  for (int i=0;i<N;i++) {
+    // get pointer to src molecule
+    molecule *srcmol = ReturnIndex(src[i]);
+    //srcmol->Center.Zero();
+    srcmol->Translate(&srcmol->Center);
+  }
+  // perform a simple multi merge
+  SimpleMultiMerge(mol, src, N);
+  return true;
+};
+
+/** Embedding merge of a given set of molecules into one.
+ * Embedding merge inserts one molecule into the other.
+ * \param *mol destination molecule
+ * \param *srcmol source molecule
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ * \TODO find embedding center
+ */
+bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+    return false;
+
+  // calculate center for merge
+  srcmol->Center.CopyVector(mol->FindEmbeddingHole((ofstream *)&cout, srcmol));
+  srcmol->Center.Zero();
+
+  // perform simple merge
+  SimpleMerge(mol, srcmol);
+  return true;
 };
 
@@ -127 +315 @@
 void MoleculeListClass::Output(ofstream *out)
 {
-  *out<< Verbose(1) << "MoleculeList: ";
-  for (int i=0;i<NumberOfMolecules;i++)
-    *out << ListOfMolecules[i] << "\t";
+  *out << Verbose(1) << "MoleculeList: ";
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
+    *out << *ListRunner << "\t";
   *out << endl;
 };
@@ -145 +333 @@
   atom *Runner = NULL;
   double ***FitConstant = NULL, **correction = NULL;
-  int a,b;
+  int a, b;
   ofstream output;
   ifstream input;
@@ -165 +353 @@
   input.open(line.c_str());
   if (input == NULL) {
-    cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+    cerr << endl << "Unable to open " << line << ", is the directory correct?"
+        << endl;
     return false;
   }
-  a=0;
-  b=-1; // we overcount by one 
+  a = 0;
+  b = -1; // we overcount by one
   while (!input.eof()) {
     input.getline(ParsedLine, 1023);
     zeile.str(ParsedLine);
-    int i=0;
+    int i = 0;
     while (!zeile.eof()) {
       zeile >> distance;
       i++;
     }
-    if (i > a) 
-      a = i; 
+    if (i > a)
+      a = i;
     b++;
   }
   cout << "I recognized " << a << " columns and " << b << " rows, ";
   input.close();
-  
+
   // 0b. allocate memory for constants
-  FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
-  for (int k=0;k<3;k++) {
-    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
-    for (int i=a;i--;) {
-      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+  FitConstant = (double ***) Malloc(sizeof(double **) * 3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k = 0; k < 3; k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i = a; i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     }
   }
   // 0c. parse in constants
-  for (int i=0;i<3;i++) {
+  for (int i = 0; i < 3; i++) {
     line = path;
     line.append("/");
     line += FRAGMENTPREFIX;
-    sprintf(ParsedLine, "%d", i+1);
+    sprintf(ParsedLine, "%d", i + 1);
     line += ParsedLine;
     line += FITCONSTANTSUFFIX;
@@ -206 +395 @@
       return false;
     }
-    int k = 0,l;
+    int k = 0, l;
     while ((!input.eof()) && (k < b)) {
       input.getline(ParsedLine, 1023);
@@ -223 +412 @@
     input.close();
   }
-  for(int k=0;k<3;k++) {
+  for (int k = 0; k < 3; k++) {
     cout << "Constants " << k << ":" << endl;
-    for (int j=0;j<b;j++) {
-      for (int i=0;i<a;i++) {
+    for (int j = 0; j < b; j++) {
+      for (int i = 0; i < a; i++) {
         cout << FitConstant[k][i][j] << "\t";
       }
@@ -233 +422 @@
     cout << endl;
   }
-  
+
   // 0d. allocate final correction matrix
-  correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
-  for (int i=a;i--;)
-    correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
-        
+  correction = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **correction");
+  for (int i = a; i--;)
+    correction[i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
+
   // 1a. go through every molecule in the list
-  for(int i=NumberOfMolecules;i--;) {
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     // 1b. zero final correction matrix
-    for (int k=a;k--;)
-      for (int j=b;j--;)
+    for (int k = a; k--;)
+      for (int j = b; j--;)
         correction[k][j] = 0.;
     // 2. take every hydrogen that is a saturated one
-    Walker = ListOfMolecules[i]->start;
-    while (Walker->next != ListOfMolecules[i]->end) {
+    Walker = (*ListRunner)->start;
+    while (Walker->next != (*ListRunner)->end) {
       Walker = Walker->next;
-      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
-      if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
-        Runner = ListOfMolecules[i]->start;
-        while (Runner->next != ListOfMolecules[i]->end) {
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      if ((Walker->type->Z == 1) && ((Walker->father == NULL)
+          || (Walker->father->type->Z != 1))) { // if it's a hydrogen
+        Runner = (*ListRunner)->start;
+        while (Runner->next != (*ListRunner)->end) {
           Runner = Runner->next;
-          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << 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) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {  // (hydrogens have only one 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!)
             // 4. evaluate the morse potential for each matrix component and add up
-            distance = sqrt(Runner->x.Distance(&Walker->x));
-            //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
-            for(int k=0;k<a;k++) {
-              for (int j=0;j<b;j++) {
-                switch(k) {
-                case 1:
-                case 7:
-                case 11:
-                  tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
-                  break;
-                default:
-                  tmp = FitConstant[0][k][j] * pow( distance,  FitConstant[1][k][j]) + FitConstant[2][k][j];
+            distance = Runner->x.Distance(&Walker->x);
+            //cout << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
+            for (int k = 0; k < a; k++) {
+              for (int j = 0; j < b; j++) {
+                switch (k) {
+                  case 1:
+                  case 7:
+                  case 11:
+                    tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
+                    break;
+                  default:
+                    tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
                 };
-                correction[k][j] -= tmp;    // ground state is actually lower (disturbed by additional interaction)
+                correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
                 //cout << tmp << "\t";
               }
@@ -285 +475 @@
     line.append("/");
     line += FRAGMENTPREFIX;
-    FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
+    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
     line += FragmentNumber;
-    delete(FragmentNumber);
+    delete (FragmentNumber);
     line += HCORRECTIONSUFFIX;
     output.open(line.c_str());
     output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
-    for (int j=0;j<b;j++) {
-      for(int i=0;i<a;i++)
+    for (int j = 0; j < b; j++) {
+      for (int i = 0; i < a; i++)
         output << correction[i][j] << "\t";
       output << endl;
@@ -303 +493 @@
   output.open(line.c_str());
   output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
-  for (int j=0;j<b;j++) {
-    for(int i=0;i<a;i++)
+  for (int j = 0; j < b; j++) {
+    for (int i = 0; i < a; i++)
       output << 0 << "\t";
     output << endl;
@@ -310 +500 @@
   output.close();
   // 6. free memory of parsed matrices
-  FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
-  for (int k=0;k<3;k++) {
-    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
-    for (int i=a;i--;) {
-      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+  FitConstant = (double ***) Malloc(sizeof(double **) * a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k = 0; k < 3; k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i = a; i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     }
   }
@@ -327 +517 @@
  * \return true - file written successfully, false - writing failed
  */
-bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex)
+bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path,
+    int *SortIndex)
 {
   bool status = true;
@@ -342 +533 @@
     //cout << Verbose(1) << "Final AtomicForcesList: ";
     //output << prefix << "Forces" << endl;
-    for(int j=0;j<NumberOfMolecules;j++) {
-      //if (TEList[j] != 0) {
-      runner = ListOfMolecules[j]->elemente->start;
-      while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
+    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+      runner = (*ListRunner)->elemente->start;
+      while (runner->next != (*ListRunner)->elemente->end) { // go through every element
         runner = runner->next;
-        if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
-          Walker = ListOfMolecules[j]->start;
-          while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
+        if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms
+          Walker = (*ListRunner)->start;
+          while (Walker->next != (*ListRunner)->end) { // go through every atom of this element
             Walker = Walker->next;
             if (Walker->type->Z == runner->Z) {
               if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
                 //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
-                ForcesFile <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
-                } else  // otherwise a -1 to indicate an added saturation hydrogen
-                  ForcesFile << "-1\t";
-              }
+                ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t";
+              } else
+                // otherwise a -1 to indicate an added saturation hydrogen
+                ForcesFile << "-1\t";
             }
           }
+        }
       }
       ForcesFile << endl;
@@ -370 +561 @@
   }
   ForcesFile.close();
-  
+
   return status;
 };
@@ -380 +571 @@
  * \return true - success (each file was written), false - something went wrong.
  */
-bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
+bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out,
+    config *configuration, int *SortIndex)
 {
   ofstream outputFragment;
@@ -393 +585 @@
   int FragmentCounter = 0;
   ofstream output;
-  string basis("3-21G");
-  
+
   // store the fragments as config and as xyz
-  for(int i=0;i<NumberOfMolecules;i++) {
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     // save default path as it is changed for each fragment
     path = configuration->GetDefaultPath();
@@ -405 +596 @@
 
     // correct periodic
-    //ListOfMolecules[i]->ScanForPeriodicCorrection(out);
+    (*ListRunner)->ScanForPeriodicCorrection(out);
 
     // output xyz file
-    FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
+    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
     sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
-    if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ...";
+    if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
       *out << " done." << endl;
     else
@@ -422 +613 @@
     // list atoms in fragment for debugging
     *out << Verbose(2) << "Contained atoms: ";
-    Walker = ListOfMolecules[i]->start;
-    while (Walker->next != ListOfMolecules[i]->end) {
+    Walker = (*ListRunner)->start;
+    while (Walker->next != (*ListRunner)->end) {
       Walker = Walker->next;
       *out << Walker->Name << " ";
     }
     *out << endl;
-    
+
     // center on edge
-    ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
-    ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
+    (*ListRunner)->CenterEdge(out, &BoxDimension);
+    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
     int j = -1;
-    for (int k=0;k<NDIM;k++) {
-      j += k+1;
-      BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
-      ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
-    }
-    ListOfMolecules[i]->Translate(&BoxDimension);
+    for (int k = 0; k < NDIM; k++) {
+      j += k + 1;
+      BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
+      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
+    }
+    (*ListRunner)->Translate(&BoxDimension);
 
     // also calculate necessary orbitals
-    ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
-    ListOfMolecules[i]->CalculateOrbitals(*configuration);
-    
+    (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
+    (*ListRunner)->CalculateOrbitals(*configuration);
+
     // change path in config
     //strcpy(PathBackup, configuration->configpath);
     sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     configuration->SetDefaultPath(FragmentName);
-    
+
     // and save as config
     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
-    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
-    if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...";
+    if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner))))
       *out << " done." << endl;
     else
@@ -461 +652 @@
     configuration->SetDefaultPath(PathBackup);
 
-
     // and save as mpqc input file
     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
-    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
-    if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...";
+    if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner))))
       *out << " done." << endl;
     else
       *out << " failed." << endl;
-       
+
     result = result && intermediateResult;
     //outputFragment.close();
     //outputFragment.clear();
-    delete(FragmentNumber);
+    delete (FragmentNumber);
     //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
   }
   cout << " done." << endl;
-  
+
   // printing final number
   *out << "Final number of fragments: " << FragmentCounter << "." << endl;
-      
+
   return result;
 };
+
+/** Counts the number of molecules with the molecule::ActiveFlag set.
+ * \return number of molecules with ActiveFlag set to true.
+ */
+int MoleculeListClass::NumberOfActiveMolecules()
+{
+  int count = 0;
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
+    count += ((*ListRunner)->ActiveFlag ? 1 : 0);
+  return count;
+};
+
 
 /******************************************* Class MoleculeLeafClass ************************************************/
@@ -493 +695 @@
 MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
 {
-//  if (Up != NULL)
-//    if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
-//      Up->DownLeaf = this;
-//  UpLeaf = Up;
-//  DownLeaf = NULL;
+  //  if (Up != NULL)
+  //    if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
+  //      Up->DownLeaf = this;
+  //  UpLeaf = Up;
+  //  DownLeaf = NULL;
   Leaf = NULL;
   previous = PreviousLeaf;
@@ -513 +715 @@
 MoleculeLeafClass::~MoleculeLeafClass()
 {
-//  if (DownLeaf != NULL) {// drop leaves further down
-//    MoleculeLeafClass *Walker = DownLeaf;
-//    MoleculeLeafClass *Next;
-//    do {
-//      Next = Walker->NextLeaf;
-//      delete(Walker);
-//      Walker = Next;
-//    } while (Walker != NULL);
-//    // Last Walker sets DownLeaf automatically to NULL
-//  }
+  //  if (DownLeaf != NULL) {// drop leaves further down
+  //    MoleculeLeafClass *Walker = DownLeaf;
+  //    MoleculeLeafClass *Next;
+  //    do {
+  //      Next = Walker->NextLeaf;
+  //      delete(Walker);
+  //      Walker = Next;
+  //    } while (Walker != NULL);
+  //    // Last Walker sets DownLeaf automatically to NULL
+  //  }
   // remove the leaf itself
   if (Leaf != NULL) {
-    delete(Leaf);
+    delete (Leaf);
     Leaf = NULL;
   }
   // remove this Leaf from level list
-  if (previous != NULL)   
+  if (previous != NULL)
     previous->next = next;
-//  } else { // we are first in list (connects to UpLeaf->DownLeaf)
-//    if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
-//      NextLeaf->UpLeaf = UpLeaf;  // either null as we are top level or the upleaf of the first node
-//    if (UpLeaf != NULL)
-//      UpLeaf->DownLeaf = NextLeaf;  // either null as we are only leaf or NextLeaf if we are just the first
-//  }
-//  UpLeaf = NULL;
+  //  } else { // we are first in list (connects to UpLeaf->DownLeaf)
+  //    if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
+  //      NextLeaf->UpLeaf = UpLeaf;  // either null as we are top level or the upleaf of the first node
+  //    if (UpLeaf != NULL)
+  //      UpLeaf->DownLeaf = NextLeaf;  // either null as we are only leaf or NextLeaf if we are just the first
+  //  }
+  //  UpLeaf = NULL;
   if (next != NULL) // are we last in list
     next->previous = previous;
@@ -576 +778 @@
     return false;
   }
-  
+
   if (status) {
-    *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
+    *out << Verbose(1) << "Creating adjacency list for subgraph " << this
+        << "." << endl;
     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
+      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
+        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->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, Binder->BondDegree);
         } else {
           *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
@@ -601 +804 @@
     FragmentCounter--;
   }
-  
+
   if ((FreeList) && (ListOfLocalAtoms != NULL)) {
     // free the index lookup list
-    Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
+    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
     if (FragmentCounter == 0) // first fragments frees the initial pointer to list
-      Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
+      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
   }
   FragmentCounter--;
@@ -621 +824 @@
  * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
  */
-bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
+bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out,
+    KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
 {
   atom *Walker = NULL, *Father = NULL;
 
   if (RootStack != NULL) {
-    // find first root candidates 
+    // find first root candidates
     if (&(RootStack[FragmentCounter]) != NULL) {
-      RootStack[FragmentCounter].clear();  
+      RootStack[FragmentCounter].clear();
       Walker = Leaf->start;
       while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
@@ -634 +838 @@
         Father = Walker->GetTrueFather();
         if (AtomMask[Father->nr]) // apply mask
-    #ifdef ADDHYDROGEN
+#ifdef ADDHYDROGEN
           if (Walker->type->Z != 1) // skip hydrogen
-    #endif
-            RootStack[FragmentCounter].push_front(Walker->nr);
+#endif
+          RootStack[FragmentCounter].push_front(Walker->nr);
       }
       if (next != NULL)
         next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
-    }  else {
-      *out << Verbose(1) << "Rootstack[" << FragmentCounter  << "] is NULL." << endl;
+    } else {
+      *out << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl;
       return false;
     }
@@ -664 +868 @@
 {
   bool status = true;
-  
+
   int Counter = Count();
   if (ListOfLocalAtoms == NULL) { // allocated initial pointer
     // allocate and set each field to NULL
-    ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
+    ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **) * Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
     if (ListOfLocalAtoms != NULL) {
-      for (int i=Counter;i--;)
+      for (int i = Counter; i--;)
         ListOfLocalAtoms[i] = NULL;
       FreeList = FreeList && true;
@@ -676 +880 @@
       status = false;
   }
-  
+
   if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
     status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
     FreeList = FreeList && true;
   }
-  
+
   return status;
 };
@@ -695 +899 @@
  * \retuen true - success, false - failure
  */
-bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
+bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out,
+    molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms,
+    Graph **&FragmentList, int &FragmentCounter, bool FreeList)
 {
   bool status = true;
   int KeySetCounter = 0;
-  
+
   *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
   // fill ListOfLocalAtoms if NULL was given
@@ -710 +916 @@
   if (FragmentList == NULL) {
     KeySetCounter = Count();
-    FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
-    for(int i=KeySetCounter;i--;)
+    FragmentList = (Graph **) Malloc(sizeof(Graph *) * KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
+    for (int i = KeySetCounter; i--;)
       FragmentList[i] = NULL;
     KeySetCounter = 0;
   }
-  
+
   if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
     // assign scanned keysets
@@ -721 +927 @@
       FragmentList[FragmentCounter] = new Graph;
     KeySet *TempSet = new KeySet;
-    for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
-      if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
+    for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers!
+      if (ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
         // translate keyset to local numbers
-        for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+        for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
           TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
         // insert into FragmentList
-        FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
+        FragmentList[FragmentCounter]->insert(GraphPair(*TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
       }
       TempSet->clear();
     }
-    delete(TempSet);
+    delete (TempSet);
     if (KeySetCounter == 0) {// if there are no keysets, delete the list
       *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
-      delete(FragmentList[FragmentCounter]);
+      delete (FragmentList[FragmentCounter]);
     } else
       *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
@@ -743 +949 @@
   } else
     *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
-  
+
   if ((FreeList) && (ListOfLocalAtoms != NULL)) {
     // free the index lookup list
-    Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
+    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
     if (FragmentCounter == 0) // first fragments frees the initial pointer to list
-      Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
+      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
   }
   *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
@@ -760 +966 @@
  * \param &TotalNumberOfKeySets global key set counter
  * \param &TotalGraph Graph to be filled with global numbers
- */ 
-void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
+ */
+void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out,
+    Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets,
+    Graph &TotalGraph)
 {
   *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
   KeySet *TempSet = new KeySet;
   if (FragmentList[FragmentCounter] != NULL) {
-    for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
-      for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+    for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
+      for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
         TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
-      TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
+      TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
       TempSet->clear();
     }
-    delete(TempSet);
+    delete (TempSet);
   } else {
     *out << Verbose(1) << "FragmentList is NULL." << endl;
@@ -788 +996 @@
 {
   if (next != NULL)
-    return next->Count()+1;
+    return next->Count() + 1;
   else
-    return 1; 
-};
+    return 1;
+};
+