Changeset ce5ac3 for src/molecules.cpp

Timestamp:

Jul 23, 2009, 12:32:48 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

d067d45

Parents:

986c80

Message:

Fix indentation from tabs to two spaces to prepare merging with MultipleMolecules

File:

• src/molecules.cpp (modified) (32 diffs)

src/molecules.cpp

@@ -53 +53 @@
   BondCount = 0;
   NoNonBonds = 0;
-        NoNonHydrogen = 0;
+  NoNonHydrogen = 0;
   NoCyclicBonds = 0;
   ListOfBondsPerAtom = NULL;
@@ -117 +117 @@
 {
   if (pointer != NULL) {
-        atom *walker = new atom();
-        walker->type = pointer->type;   // copy element of atom
-        walker->x.CopyVector(&pointer->x); // copy coordination
+    atom *walker = new atom();
+    walker->type = pointer->type;  // copy element of atom
+    walker->x.CopyVector(&pointer->x); // copy coordination
     walker->v.CopyVector(&pointer->v); // copy velocity
     walker->FixedIon = pointer->FixedIon;
@@ -186 +186 @@
   double *matrix;
 
-//      *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
+//  *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
   InBondvector.CopyVector(&TopReplacement->x);
@@ -427 +427 @@
   }
 
-//      *out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
+//  *out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
   return AllWentWell;
 };
@@ -753 +753 @@
 Vector * molecule::DetermineCenterOfGravity(ofstream *out)
 {
-        atom *ptr = start->next;  // start at first in list
-        Vector *a = new Vector();
-        Vector tmp;
+  atom *ptr = start->next;  // start at first in list
+  Vector *a = new Vector();
+  Vector tmp;
   double Num = 0;
 
-        a->Zero();
+  a->Zero();
 
   if (ptr != end) {   //list not empty?
@@ -909 +909 @@
 void molecule::PrincipalAxisSystem(ofstream *out, bool DoRotate)
 {
-        atom *ptr = start;  // start at first in list
-        double InertiaTensor[NDIM*NDIM];
-        Vector *CenterOfGravity = DetermineCenterOfGravity(out);
-
-        CenterGravity(out, CenterOfGravity);
-
-        // reset inertia tensor
-        for(int i=0;i<NDIM*NDIM;i++)
-                InertiaTensor[i] = 0.;
-
-        // sum up inertia tensor
-        while (ptr->next != end) {
-                ptr = ptr->next;
-                Vector x;
-                x.CopyVector(&ptr->x);
-                //x.SubtractVector(CenterOfGravity);
-                InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-                InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-                InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-                InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-                InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-                InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-                InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-                InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-                InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
-        }
-        // print InertiaTensor for debugging
-        *out << "The inertia tensor is:" << endl;
-        for(int i=0;i<NDIM;i++) {
-          for(int j=0;j<NDIM;j++)
-            *out << InertiaTensor[i*NDIM+j] << " ";
-          *out << endl;
-        }
-        *out << endl;
-
-        // diagonalize to determine principal axis system
-        gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
-        gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
-        gsl_vector *eval = gsl_vector_alloc(NDIM);
-        gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
-        gsl_eigen_symmv(&m.matrix, eval, evec, T);
-        gsl_eigen_symmv_free(T);
-        gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
-
-        for(int i=0;i<NDIM;i++) {
-                *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
-                *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
-        }
-
-        // check whether we rotate or not
-        if (DoRotate) {
-          *out << Verbose(1) << "Transforming molecule into PAS ... ";
-          // the eigenvectors specify the transformation matrix
-          ptr = start;
-          while (ptr->next != end) {
-            ptr = ptr->next;
-            for (int j=0;j<MDSteps;j++)
-              Trajectories[ptr].R.at(j).MatrixMultiplication(evec->data);
-            ptr->x.MatrixMultiplication(evec->data);
-          }
-          *out << "done." << endl;
-
-          // summing anew for debugging (resulting matrix has to be diagonal!)
-          // reset inertia tensor
+  atom *ptr = start;  // start at first in list
+  double InertiaTensor[NDIM*NDIM];
+  Vector *CenterOfGravity = DetermineCenterOfGravity(out);
+
+  CenterGravity(out, CenterOfGravity);
+
+  // reset inertia tensor
+  for(int i=0;i<NDIM*NDIM;i++)
+    InertiaTensor[i] = 0.;
+
+  // sum up inertia tensor
+  while (ptr->next != end) {
+    ptr = ptr->next;
+    Vector x;
+    x.CopyVector(&ptr->x);
+    //x.SubtractVector(CenterOfGravity);
+    InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
+    InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
+    InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
+    InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
+    InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
+    InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
+    InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
+    InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
+    InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+  }
+  // print InertiaTensor for debugging
+  *out << "The inertia tensor is:" << endl;
+  for(int i=0;i<NDIM;i++) {
+    for(int j=0;j<NDIM;j++)
+      *out << InertiaTensor[i*NDIM+j] << " ";
+    *out << endl;
+  }
+  *out << endl;
+
+  // diagonalize to determine principal axis system
+  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
+  gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
+  gsl_vector *eval = gsl_vector_alloc(NDIM);
+  gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
+  gsl_eigen_symmv(&m.matrix, eval, evec, T);
+  gsl_eigen_symmv_free(T);
+  gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
+
+  for(int i=0;i<NDIM;i++) {
+    *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
+    *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
+  }
+
+  // check whether we rotate or not
+  if (DoRotate) {
+    *out << Verbose(1) << "Transforming molecule into PAS ... ";
+    // the eigenvectors specify the transformation matrix
+    ptr = start;
+    while (ptr->next != end) {
+      ptr = ptr->next;
+      for (int j=0;j<MDSteps;j++)
+        Trajectories[ptr].R.at(j).MatrixMultiplication(evec->data);
+      ptr->x.MatrixMultiplication(evec->data);
+    }
+    *out << "done." << endl;
+
+    // summing anew for debugging (resulting matrix has to be diagonal!)
+    // reset inertia tensor
     for(int i=0;i<NDIM*NDIM;i++)
       InertiaTensor[i] = 0.;
@@ -1001 +1001 @@
     }
     *out << endl;
-        }
-
-        // free everything
-        delete(CenterOfGravity);
-        gsl_vector_free(eval);
-        gsl_matrix_free(evec);
+  }
+
+  // free everything
+  delete(CenterOfGravity);
+  gsl_vector_free(eval);
+  gsl_matrix_free(evec);
 };
 
@@ -2027 +2027 @@
     runner = elemente->start;
     while (runner->next != elemente->end) { // go through every element
-                runner = runner->next;
+      runner = runner->next;
       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
         ElementNo++;
@@ -2121 +2121 @@
 bool molecule::Checkout(ofstream *out)  const
 {
-        return elemente->Checkout(out, ElementsInMolecule);
+  return elemente->Checkout(out, ElementsInMolecule);
 };
 
@@ -2174 +2174 @@
     runner = elemente->start;
     while (runner->next != elemente->end) { // go through every element
-                runner = runner->next;
+      runner = runner->next;
       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
         ElementNo++;
@@ -2196 +2196 @@
 void molecule::CountAtoms(ofstream *out)
 {
-        int i = 0;
-        atom *Walker = start;
+  int i = 0;
+  atom *Walker = start;
   while (Walker->next != end) {
     Walker = Walker->next;
@@ -2204 +2204 @@
   if ((AtomCount == 0) || (i != AtomCount)) {
     *out << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
-        AtomCount = i;
+    AtomCount = i;
 
     // count NonHydrogen atoms and give each atom a unique name
     if (AtomCount != 0) {
-            i=0;
-            NoNonHydrogen = 0;
+      i=0;
+      NoNonHydrogen = 0;
       Walker = start;
-            while (Walker->next != end) {
-              Walker = Walker->next;
-              Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-              if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
-                NoNonHydrogen++;
-              Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
-              Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
-              sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
+      while (Walker->next != end) {
+        Walker = Walker->next;
+        Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
+        if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
+          NoNonHydrogen++;
+        Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
+        Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
+        sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
         *out << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
-              i++;
-            }
+        i++;
+      }
     } else
-        *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
+      *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
   }
 };
@@ -2231 +2231 @@
 void molecule::CountElements()
 {
-        int i = 0;
+  int i = 0;
   for(i=MAX_ELEMENTS;i--;)
-        ElementsInMolecule[i] = 0;
-        ElementCount = 0;
+    ElementsInMolecule[i] = 0;
+  ElementCount = 0;
 
   atom *walker = start;
@@ -2243 +2243 @@
   }
   for(i=MAX_ELEMENTS;i--;)
-        ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
+    ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
 };
 
@@ -2334 +2334 @@
 {
 
-        // 1 We will parse bonds out of the dbond file created by tremolo.
-                        int atom1, atom2, temp;
-                        atom *Walker, *OtherWalker;
-
-                if (!input)
-                {
-                        cout << Verbose(1) << "Opening silica failed \n";
-                };
-
-                        *input >> ws >> atom1;
-                        *input >> ws >> atom2;
-                cout << Verbose(1) << "Scanning file\n";
-                while (!input->eof()) // Check whether we read everything already
-                {
-                                *input >> ws >> atom1;
-                                *input >> ws >> atom2;
-                        if(atom2<atom1) //Sort indices of atoms in order
-                        {
-                                temp=atom1;
-                                atom1=atom2;
-                                atom2=temp;
-                        };
-
-                        Walker=start;
-                        while(Walker-> nr != atom1) // Find atom corresponding to first index
-                        {
-                                Walker = Walker->next;
-                        };
-                        OtherWalker = Walker->next;
-                        while(OtherWalker->nr != atom2) // Find atom corresponding to second index
-                        {
-                                OtherWalker= OtherWalker->next;
-                        };
-                        AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
-                        
-                }
-
-                CreateListOfBondsPerAtom(out);
+  // 1 We will parse bonds out of the dbond file created by tremolo.
+      int atom1, atom2, temp;
+      atom *Walker, *OtherWalker;
+
+          if (!input)
+          {
+            cout << Verbose(1) << "Opening silica failed \n";
+          };
+
+      *input >> ws >> atom1;
+      *input >> ws >> atom2;
+          cout << Verbose(1) << "Scanning file\n";
+          while (!input->eof()) // Check whether we read everything already
+          {
+        *input >> ws >> atom1;
+        *input >> ws >> atom2;
+            if(atom2<atom1) //Sort indices of atoms in order
+            {
+              temp=atom1;
+              atom1=atom2;
+              atom2=temp;
+            };
+
+            Walker=start;
+            while(Walker-> nr != atom1) // Find atom corresponding to first index
+            {
+              Walker = Walker->next;
+            };
+            OtherWalker = Walker->next;
+            while(OtherWalker->nr != atom2) // Find atom corresponding to second index
+            {
+              OtherWalker= OtherWalker->next;
+            };
+            AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
+            
+          }
+
+          CreateListOfBondsPerAtom(out);
 
 };
@@ -2520 +2520 @@
     // preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
     // double bonds as was expected.
-                if (BondCount != 0) {
+    if (BondCount != 0) {
       NoCyclicBonds = 0;
-                  *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
-                  do {
-                    No = 0; // No acts as breakup flag (if 1 we still continue)
+      *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
+      do {
+        No = 0; // No acts as breakup flag (if 1 we still continue)
         Walker = start;
         while (Walker->next != end) { // go through every atom
@@ -2538 +2538 @@
             for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
               OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-                    // count valence of second partner
+              // count valence of second partner
               NoBonds = 0;
               for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
                 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
               *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-                    if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
-                      if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
+              if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
+                if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
                   Candidate = OtherWalker;
                   CandidateBondNo = i;
-                        *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
-                      }
-                    }
-                  }
+                  *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
+                }
+              }
+            }
             if ((Candidate != NULL) && (CandidateBondNo != -1)) {
               ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
@@ -2558 +2558 @@
               FalseBondDegree++;
           }
-                    }
-                  } while (No);
-                *out << " done." << endl;
-                } else
-                        *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
-          *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
-                
-          // output bonds for debugging (if bond chain list was correctly installed)
-          *out << Verbose(1) << endl << "From contents of bond chain list:";
-          bond *Binder = first;
+        }
+      } while (No);
+    *out << " done." << endl;
+    } else
+      *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
+    *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
+    
+    // output bonds for debugging (if bond chain list was correctly installed)
+    *out << Verbose(1) << endl << "From contents of bond chain list:";
+    bond *Binder = first;
     while(Binder->next != last) {
       Binder = Binder->next;
-                        *out << *Binder << "\t" << endl;
+      *out << *Binder << "\t" << endl;
     }
     *out << endl;
   } else
-        *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
+    *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
   Free((void **)&matrix, "molecule::CreateAdjacencyList: *matrix");
@@ -3600 +3600 @@
 int molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
 {
-        MoleculeListClass *BondFragments = NULL;
+  MoleculeListClass *BondFragments = NULL;
   int *SortIndex = NULL;
   int *MinimumRingSize = new int[AtomCount];
@@ -3806 +3806 @@
     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
     if (Walker != NULL) // if this Walker exists in the subgraph ...
-        for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
-        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-              if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
-                LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
+      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
+        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+        if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
+            LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
                 *out << Verbose(3) << "Found local edge " << *(ListOfBondsPerAtom[Walker->nr][i]) << "." << endl;
-                break;
-            }
-        }
+          break;
+        }
+      }
     Binder = ReferenceStack->PopFirst();  // loop the stack for next item
     *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
@@ -4308 +4308 @@
   for (int i=0;i<AtomCount;i++) { // reset all atom labels
     // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
-                Labels[i] = -1;
+    Labels[i] = -1;
     SonList[i] = NULL;
     PredecessorList[i] = NULL;
@@ -4316 +4316 @@
   for (int i=0;i<BondCount;i++)
     ColorEdgeList[i] = white;
-        RootStack->ClearStack();        // clearstack and push first atom if exists
+  RootStack->ClearStack();  // clearstack and push first atom if exists
   TouchedStack->ClearStack();
   Walker = start->next;
@@ -4335 +4335 @@
   ///// OUTER LOOP ////////////
   while (!RootStack->IsEmpty()) {
-        // get new root vertex from atom stack
-        Root = RootStack->PopFirst();
+    // get new root vertex from atom stack
+    Root = RootStack->PopFirst();
     ShortestPathList[Root->nr] = 0;
     if (Labels[Root->nr] == -1)
@@ -4344 +4344 @@
     *out << Verbose(0) << "Root for this loop is: " << Root->Name << ".\n";
 
-                // clear snake stack
-          SnakeStack->ClearStack();
+    // clear snake stack
+    SnakeStack->ClearStack();
     //SnakeStack->TestImplementation(out, start->next);
 
@@ -4352 +4352 @@
     // - what about cyclic bonds?
     Walker = Root;
-        do {
+    do {
       *out << Verbose(1) << "Current Walker is: " << Walker->Name;
       // initial setting of the new Walker: label, color, shortest path and put on stacks
-                if (Labels[Walker->nr] == -1)   {       // give atom a unique, monotonely increasing number
-                        Labels[Walker->nr] = RunningIndex++;
+      if (Labels[Walker->nr] == -1)  {  // give atom a unique, monotonely increasing number
+        Labels[Walker->nr] = RunningIndex++;
         RootStack->Push(Walker);
       }
       *out << ", has label " << Labels[Walker->nr];
-                if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) {     // color it if newly discovered and push on stacks (and if within reach!)
+      if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) {  // color it if newly discovered and push on stacks (and if within reach!)
         if ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white)) {
           // Binder ought to be set still from last neighbour search
@@ -4374 +4374 @@
           ColorVertexList[Walker->nr] = darkgray; // mark as dark gray of on snake stack
         }
-                }
+      }
       *out << ", SP of " << ShortestPathList[Walker->nr]  << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
 
       // then check the stack for a newly stumbled upon fragment
-                if (SnakeStack->ItemCount() == Order) { // is stack full?
+      if (SnakeStack->ItemCount() == Order) { // is stack full?
         // store the fragment if it is one and get a removal candidate
         Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
@@ -4391 +4391 @@
           }
         }
-                } else
+      } else
         Removal = NULL;
 
@@ -4398 +4398 @@
       if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {  // don't look, if a new walker has been set above
         *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
-                        OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
+        OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
         if (ShortestPathList[Walker->nr] < Order) {
-                        for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
+          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
             Binder = ListOfBondsPerAtom[Walker->nr][i];
             *out << Verbose(2) << "Current bond is " << *Binder << ": ";
-                                OtherAtom = Binder->GetOtherAtom(Walker);
+            OtherAtom = Binder->GetOtherAtom(Walker);
             if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
               *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
               //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
             } else { // otherwise check its colour and element
-                                if (
+              if (
 #ifdef ADDHYDROGEN
               (OtherAtom->type->Z != 1) &&
@@ -4430 +4430 @@
               }
             }
-                        }
+          }
         } else {  // means we have stepped beyond the horizon: Return!
           Walker = PredecessorList[Walker->nr];
@@ -4436 +4436 @@
           *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
         }
-                        if (Walker != OtherAtom) {      // if no white neighbours anymore, color it black
+        if (Walker != OtherAtom) {  // if no white neighbours anymore, color it black
           *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
-                                ColorVertexList[Walker->nr] = black;
-                                Walker = PredecessorList[Walker->nr];
-                        }
-      }
-        } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
+          ColorVertexList[Walker->nr] = black;
+          Walker = PredecessorList[Walker->nr];
+        }
+      }
+    } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
     *out << Verbose(2) << "Inner Looping is finished." << endl;
 
@@ -4562 +4562 @@
         bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
         if (bit) {  // if bit is set, we add this bond partner
-                OtherWalker = BondsSet[j]->rightatom;   // rightatom is always the one more distant, i.e. the one to add
+          OtherWalker = BondsSet[j]->rightatom;  // rightatom is always the one more distant, i.e. the one to add
           //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
@@ -4584 +4584 @@
         if (SubOrder > 1) {    // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
           // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
-          SP = RootDistance+1;  // this is the next level
+          SP = RootDistance+1;  // this is the next level
           // first count the members in the subset
           SubSetDimension = 0;
-          Binder = FragmentSearch->BondsPerSPList[2*SP];                // start node for this level
-          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {              // compare to end node of this level
+          Binder = FragmentSearch->BondsPerSPList[2*SP];    // start node for this level
+          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {    // compare to end node of this level
             Binder = Binder->next;
             for (int k=TouchedIndex;k--;) {