Changes in molecuilder/src/moleculelist.cpp [1b2aa1:f23d7d]

File:

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

molecuilder/src/moleculelist.cpp

@@ -1 +1 @@
 /** \file MoleculeListClass.cpp
- *
+ * 
  * Function implementations for the class MoleculeListClass.
- *
+ * 
  */
 
@@ -13 +13 @@
 MoleculeListClass::MoleculeListClass()
 {
-  // empty lists
-  ListOfMolecules.clear();
-  MaxIndex = 1;
-};
+};
+
+/** 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;
+};
+
 
 /** Destructor for MoleculeListClass.
@@ -23 +34 @@
 {
   cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
-  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
-    cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl;
-    delete (*ListRunner);
+  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;
+    }
   }
   cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
-  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);
+  Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
 };
 
@@ -53 +57 @@
   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;
@@ -87 +90 @@
       // 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;
@@ -111 +114 @@
         }
       }
-      delete[] (aList);
-      delete[] (bList);
+      delete[](aList);
+      delete[](bList);
       return flag;
     }
   }
-  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;
+  return  -1;
 };
 
@@ -315 +127 @@
 void MoleculeListClass::Output(ofstream *out)
 {
-  *out << Verbose(1) << "MoleculeList: ";
-  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
-    *out << *ListRunner << "\t";
+  *out<< Verbose(1) << "MoleculeList: ";
+  for (int i=0;i<NumberOfMolecules;i++)
+    *out << ListOfMolecules[i] << "\t";
   *out << endl;
 };
@@ -333 +145 @@
   atom *Runner = NULL;
   double ***FitConstant = NULL, **correction = NULL;
-  int a, b;
+  int a,b;
   ofstream output;
   ifstream input;
@@ -353 +165 @@
   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;
@@ -395 +206 @@
       return false;
     }
-    int k = 0, l;
+    int k = 0,l;
     while ((!input.eof()) && (k < b)) {
       input.getline(ParsedLine, 1023);
@@ -412 +223 @@
     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";
       }
@@ -422 +233 @@
     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 (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+  for(int i=NumberOfMolecules;i--;) {
     // 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 = (*ListRunner)->start;
-    while (Walker->next != (*ListRunner)->end) {
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
       Walker = Walker->next;
-      //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) {
+      //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) {
           Runner = Runner->next;
-          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->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) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one 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!)
             // 4. evaluate the morse potential for each matrix component and add up
-            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];
+            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];
                 };
-                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";
               }
@@ -475 +285 @@
     line.append("/");
     line += FRAGMENTPREFIX;
-    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
     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;
@@ -493 +303 @@
   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;
@@ -500 +310 @@
   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[][]");
     }
   }
@@ -517 +327 @@
  * \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;
@@ -533 +342 @@
     //cout << Verbose(1) << "Final AtomicForcesList: ";
     //output << prefix << "Forces" << endl;
-    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
-      runner = (*ListRunner)->elemente->start;
-      while (runner->next != (*ListRunner)->elemente->end) { // go through every element
+    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
         runner = runner->next;
-        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
+        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
             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;
@@ -561 +370 @@
   }
   ForcesFile.close();
-
+  
   return status;
 };
@@ -571 +380 @@
  * \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;
@@ -585 +393 @@
   int FragmentCounter = 0;
   ofstream output;
-
+  string basis("3-21G");
+  
   // store the fragments as config and as xyz
-  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+  for(int i=0;i<NumberOfMolecules;i++) {
     // save default path as it is changed for each fragment
     path = configuration->GetDefaultPath();
@@ -596 +405 @@
 
     // correct periodic
-    (*ListRunner)->ScanForPeriodicCorrection(out);
+    //ListOfMolecules[i]->ScanForPeriodicCorrection(out);
 
     // output xyz file
-    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, 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 = (*ListRunner)->OutputXYZ(&outputFragment)))
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
+    if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
       *out << " done." << endl;
     else
@@ -613 +422 @@
     // list atoms in fragment for debugging
     *out << Verbose(2) << "Contained atoms: ";
-    Walker = (*ListRunner)->start;
-    while (Walker->next != (*ListRunner)->end) {
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
       Walker = Walker->next;
       *out << Walker->Name << " ";
     }
     *out << endl;
-
+    
     // center on edge
-    (*ListRunner)->CenterEdge(out, &BoxDimension);
-    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
+    ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
+    ListOfMolecules[i]->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);
-      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
-    }
-    (*ListRunner)->Translate(&BoxDimension);
+    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);
 
     // also calculate necessary orbitals
-    (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
-    (*ListRunner)->CalculateOrbitals(*configuration);
-
+    ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
+    ListOfMolecules[i]->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, (*ListRunner)->elemente, (*ListRunner))))
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
+    if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
       *out << " done." << endl;
     else
@@ -652 +461 @@
     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, (*ListRunner))))
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
+    if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
       *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 ************************************************/
@@ -695 +493 @@
 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;
@@ -715 +513 @@
 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;
@@ -778 +576 @@
     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;
@@ -804 +601 @@
     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--;
@@ -824 +621 @@
  * \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
@@ -838 +634 @@
         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;
     }
@@ -868 +664 @@
 {
   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;
@@ -880 +676 @@
       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;
 };
@@ -899 +695 @@
  * \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
@@ -916 +710 @@
   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
@@ -927 +721 @@
       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;
@@ -949 +743 @@
   } 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;
@@ -966 +760 @@
  * \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;
@@ -996 +788 @@
 {
   if (next != NULL)
-    return next->Count() + 1;
+    return next->Count()+1;
   else
-    return 1;
-};
-
+    return 1; 
+};