Changeset ce5ac3

Timestamp:

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

Author:

Frederik Heber <heber@…>

Branches:

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

Children:

d067d45

Parents:

986c80

Message:

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

Location:

src

Files:

• atom.cpp (modified) (2 diffs)
• bond.cpp (modified) (2 diffs)
• boundary.cpp (modified) (10 diffs)
• builder.cpp (modified) (3 diffs)
• config.cpp (modified) (1 diff)
• datacreator.cpp (modified) (1 diff)
• helpers.cpp (modified) (1 diff)
• molecules.cpp (modified) (32 diffs)
• parser.cpp (modified) (2 diffs)
• vector.cpp (modified) (6 diffs)
• verbose.cpp (modified) (3 diffs)

src/atom.cpp

@@ -13 +13 @@
 atom::atom() 
 {
-        Name = NULL;
+  Name = NULL;
   previous = NULL;
   next = NULL;
@@ -35 +35 @@
 atom::~atom() 
 {
-        Free((void **)&Name, "atom::~atom: *Name");
+  Free((void **)&Name, "atom::~atom: *Name");
   Free((void **)&ComponentNr, "atom::~atom: *ComponentNr");
 };

src/bond.cpp

@@ -14 +14 @@
 bond::bond()
 {
-        leftatom = NULL;
-        rightatom = NULL;
+  leftatom = NULL;
+  rightatom = NULL;
   previous = NULL;
   next = NULL;
-        nr = -1;
-        HydrogenBond = 0;
-        BondDegree = 0;
+  nr = -1;
+  HydrogenBond = 0;
+  BondDegree = 0;
   Used = white;
   Cyclic = false;
@@ -34 +34 @@
 bond::bond(atom *left, atom *right, int degree=1, int number=0)
 {
-        leftatom = left;
-        rightatom = right;
+  leftatom = left;
+  rightatom = right;
   previous = NULL;
   next = NULL;
-        HydrogenBond = 0;
+  HydrogenBond = 0;
   if ((left != NULL) && (right != NULL)) {
-        if ((left->type != NULL) && (left->type->Z == 1))
-                HydrogenBond++;
-        if ((right->type != NULL) && (right->type->Z == 1))
-                HydrogenBond++;
+    if ((left->type != NULL) && (left->type->Z == 1))
+      HydrogenBond++;
+    if ((right->type != NULL) && (right->type->Z == 1))
+      HydrogenBond++;
   }
   BondDegree = degree;

src/boundary.cpp

@@ -86 +86 @@
 BoundaryLineSet::~BoundaryLineSet()
 {
-        for (int i = 0; i < 2; i++) {
-                cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-                endpoints[i]->lines.erase(Nr);
-                LineMap::iterator tester = endpoints[i]->lines.begin();
-                tester++;
-                if (tester == endpoints[i]->lines.end()) {
-                        cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-                        if (endpoints[i] != NULL) {
-                                delete(endpoints[i]);
-                                endpoints[i] = NULL;
-                        } else
-                                cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
-                } else
-                        cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
-        }
+  for (int i = 0; i < 2; i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+    endpoints[i]->lines.erase(Nr);
+    LineMap::iterator tester = endpoints[i]->lines.begin();
+    tester++;
+    if (tester == endpoints[i]->lines.end()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      if (endpoints[i] != NULL) {
+        delete(endpoints[i]);
+        endpoints[i] = NULL;
+      } else
+        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
 }
 ;
@@ -181 +181 @@
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
-        for (int i = 0; i < 3; i++) {
-                cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
-                lines[i]->triangles.erase(Nr);
-                if (lines[i]->triangles.empty()) {
-                        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-                        if (lines[i] != NULL) {
-                                delete (lines[i]);
-                                lines[i] = NULL;
-                        } else
-                                cerr << "ERROR: This line " << i << " has already been free'd." << endl;
-                } else
-                        cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
-        }
+  for (int i = 0; i < 3; i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    if (lines[i]->triangles.empty()) {
+      cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+      if (lines[i] != NULL) {
+        delete (lines[i]);
+        lines[i] = NULL;
+      } else
+        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+  }
 }
 ;
@@ -573 +573 @@
 void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
 {
-        atom *Walker = mol->start;
-        bond *Binder = mol->first;
-        int i;
-        Vector *center = mol->DetermineCenterOfAll(out);
-        if (rasterfile != NULL) {
-                //cout << Verbose(1) << "Writing Raster3D file ... ";
-                *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
-                *rasterfile << "@header.r3d" << endl;
-                *rasterfile << "# All atoms as spheres" << endl;
-                while (Walker->next != mol->end) {
-                        Walker = Walker->next;
-                        *rasterfile << "2" << endl << "  ";     // 2 is sphere type
-                        for (i=0;i<NDIM;i++)
-                                *rasterfile << Walker->x.x[i]+center->x[i] << " ";
-                        *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
-                }
-
-                *rasterfile << "# All bonds as vertices" << endl;
-                while (Binder->next != mol->last) {
-                        Binder = Binder->next;
-                        *rasterfile << "3" << endl << "  ";     // 2 is round-ended cylinder type
-                        for (i=0;i<NDIM;i++)
-                                *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
-                        *rasterfile << "\t0.03\t";
-                        for (i=0;i<NDIM;i++)
-                                *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
-                        *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
-                }
-
-                *rasterfile << "# All tesselation triangles" << endl;
-                for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-                        *rasterfile << "1" << endl << "  ";     // 1 is triangle type
-                        for (i=0;i<3;i++) {     // print each node
-                                for (int j=0;j<NDIM;j++)        // and for each node all NDIM coordinates
-                                        *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
-                                *rasterfile << "\t";
-                        }
-                        *rasterfile << "1. 0. 0." << endl;      // red as colour
-                        *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
-                }
-        } else {
-                cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
-        }
-        delete(center);
+  atom *Walker = mol->start;
+  bond *Binder = mol->first;
+  int i;
+  Vector *center = mol->DetermineCenterOfAll(out);
+  if (rasterfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+    *rasterfile << "@header.r3d" << endl;
+    *rasterfile << "# All atoms as spheres" << endl;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Walker->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
+
+    *rasterfile << "# All bonds as vertices" << endl;
+    while (Binder->next != mol->last) {
+      Binder = Binder->next;
+      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t";
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
+
+    *rasterfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
+      for (i=0;i<3;i++) {  // print each node
+        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+        *rasterfile << "\t";
+      }
+      *rasterfile << "1. 0. 0." << endl;  // red as colour
+      *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  delete(center);
 };
 
@@ -946 +946 @@
 Tesselation::~Tesselation()
 {
-        cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-        for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
-                if (runner->second != NULL) {
-                        delete (runner->second);
-                        runner->second = NULL;
-                } else
-                        cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
-        }
-        for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {
-                if (runner->second != NULL) {
-                        delete (runner->second);
-                        runner->second = NULL;
-                } else
-                        cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;
-        }
-        for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-                if (runner->second != NULL) {
-                        delete (runner->second);
-                        runner->second = NULL;
-                } else
-                        cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;
-        }
+  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+  }
+  for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;
+  }
+  for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;
+  }
 }
 ;
@@ -1585 +1585 @@
     atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
 {
-        cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-        cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
-        cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
-        cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
-        cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
-        cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
-        /* OldNormal is normal vector on the old triangle
-         * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
-         * Chord points from b to a!!!
-         */
-        Vector dif_a; //Vector from a to candidate
-        Vector dif_b; //Vector from b to candidate
-        Vector AngleCheck;
-        Vector TempNormal, Umkreismittelpunkt;
-        Vector Mittelpunkt;
-
-        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
-        double BallAngle, AlternativeSign;
-        atom *Walker; // variable atom point
-
-        Vector NewUmkreismittelpunkt;
-
-        if (a != Candidate and b != Candidate and c != Candidate) {
-                cout << Verbose(3) << "We have a unique candidate!" << endl;
-                dif_a.CopyVector(&(a->x));
-                dif_a.SubtractVector(&(Candidate->x));
-                dif_b.CopyVector(&(b->x));
-                dif_b.SubtractVector(&(Candidate->x));
-                AngleCheck.CopyVector(&(Candidate->x));
-                AngleCheck.SubtractVector(&(a->x));
-                AngleCheck.ProjectOntoPlane(Chord);
-
-                SideA = dif_b.Norm();
-                SideB = dif_a.Norm();
-                SideC = Chord->Norm();
-                //Chord->Scale(-1);
-
-                alpha = Chord->Angle(&dif_a);
-                beta = M_PI - Chord->Angle(&dif_b);
-                gamma = dif_a.Angle(&dif_b);
-
-                cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-
-                if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
-                        cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-                        cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-                }
-
-                if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
-                        Umkreisradius = SideA / 2.0 / sin(alpha);
-                        //cout << Umkreisradius << endl;
-                        //cout << SideB / 2.0 / sin(beta) << endl;
-                        //cout << SideC / 2.0 / sin(gamma) << endl;
-
-                        if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
-                                cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
-                                cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
-                                sign = AngleCheck.ScalarProduct(direction1);
-                                if (fabs(sign)<MYEPSILON) {
-                                        if (AngleCheck.ScalarProduct(OldNormal)<0) {
-                                                sign =0;
-                                                AlternativeSign=1;
-                                        } else {
-                                                sign =0;
-                                                AlternativeSign=-1;
-                                        }
-                                } else {
-                                        sign /= fabs(sign);
-                                        cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
-                                }
-
-                                Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
-
-                                AngleCheck.CopyVector(&ReferencePoint);
-                                AngleCheck.Scale(-1);
-                                //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-                                AngleCheck.AddVector(&Mittelpunkt);
-                                //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-                                cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
-
-                                BallAngle = AngleCheck.Angle(OldNormal);
-                                cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-
-                                //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
-                                //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-
-                                cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
-
-                                NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
-
-                                if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
-                                        if (Storage[0]< -1.5) { // first Candidate at all
-                                                if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
-                                                        cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
-                                                        Opt_Candidate = Candidate;
-                                                        Storage[0] = sign;
-                                                        Storage[1] = AlternativeSign;
-                                                        Storage[2] = BallAngle;
-                                                        cout << " angle " << Storage[2] << " and Up/Down "
-                                                        << Storage[0] << endl;
-                                                } else
-                                                        cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
-                                        } else {
-                                                if ( Storage[2] > BallAngle) {
-                                                        if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
-                                                                cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
-                                                                Opt_Candidate = Candidate;
-                                                                Storage[0] = sign;
-                                                                Storage[1] = AlternativeSign;
-                                                                Storage[2] = BallAngle;
-                                                                cout << " angle " << Storage[2] << " and Up/Down "
-                                                                << Storage[0] << endl;
-                                                        } else
-                                                                cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
-                                                } else {
-                                                        if (DEBUG) {
-                                                                cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
-                                                        }
-                                                }
-                                        }
-                                } else {
-                                        if (DEBUG) {
-                                                cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
-                                        }
-                                }
-                        } else {
-                                if (DEBUG) {
-                                        cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
-                                }
-                        }
-                } else {
-                        if (DEBUG) {
-                                cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
-                        }
-                }
-        } else {
-                if (DEBUG) {
-                        cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
-                }
-        }
-
-        if (RecursionLevel < 5) { // Seven is the recursion level threshold.
-                for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
-                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                        if (Walker == Parent) { // don't go back the same bond
-                                continue;
-                        } else {
-                                Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
-                        }
-                }
-        }
-        cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
+  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
+  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
+  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
+  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
+  /* OldNormal is normal vector on the old triangle
+   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck;
+  Vector TempNormal, Umkreismittelpunkt;
+  Vector Mittelpunkt;
+
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
+  double BallAngle, AlternativeSign;
+  atom *Walker; // variable atom point
+
+  Vector NewUmkreismittelpunkt;
+
+  if (a != Candidate and b != Candidate and c != Candidate) {
+    cout << Verbose(3) << "We have a unique candidate!" << endl;
+    dif_a.CopyVector(&(a->x));
+    dif_a.SubtractVector(&(Candidate->x));
+    dif_b.CopyVector(&(b->x));
+    dif_b.SubtractVector(&(Candidate->x));
+    AngleCheck.CopyVector(&(Candidate->x));
+    AngleCheck.SubtractVector(&(a->x));
+    AngleCheck.ProjectOntoPlane(Chord);
+
+    SideA = dif_b.Norm();
+    SideB = dif_a.Norm();
+    SideC = Chord->Norm();
+    //Chord->Scale(-1);
+
+    alpha = Chord->Angle(&dif_a);
+    beta = M_PI - Chord->Angle(&dif_b);
+    gamma = dif_a.Angle(&dif_b);
+
+    cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+
+    if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
+      cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+      cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+    }
+
+    if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
+      Umkreisradius = SideA / 2.0 / sin(alpha);
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
+
+      if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
+        cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
+        cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
+        sign = AngleCheck.ScalarProduct(direction1);
+        if (fabs(sign)<MYEPSILON) {
+          if (AngleCheck.ScalarProduct(OldNormal)<0) {
+            sign =0;
+            AlternativeSign=1;
+          } else {
+            sign =0;
+            AlternativeSign=-1;
+          }
+        } else {
+          sign /= fabs(sign);
+          cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
+        }
+
+        Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+
+        AngleCheck.CopyVector(&ReferencePoint);
+        AngleCheck.Scale(-1);
+        //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+        AngleCheck.AddVector(&Mittelpunkt);
+        //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+        cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
+
+        BallAngle = AngleCheck.Angle(OldNormal);
+        cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+
+        //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+        //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+        cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+
+        NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
+
+        if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
+          if (Storage[0]< -1.5) { // first Candidate at all
+            if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+              cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
+              Opt_Candidate = Candidate;
+              Storage[0] = sign;
+              Storage[1] = AlternativeSign;
+              Storage[2] = BallAngle;
+              cout << " angle " << Storage[2] << " and Up/Down "
+              << Storage[0] << endl;
+            } else
+              cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+          } else {
+            if ( Storage[2] > BallAngle) {
+              if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+                cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
+                Opt_Candidate = Candidate;
+                Storage[0] = sign;
+                Storage[1] = AlternativeSign;
+                Storage[2] = BallAngle;
+                cout << " angle " << Storage[2] << " and Up/Down "
+                << Storage[0] << endl;
+              } else
+                cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+            } else {
+              if (DEBUG) {
+                cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
+              }
+            }
+          }
+        } else {
+          if (DEBUG) {
+            cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
+          }
+        }
+      } else {
+        if (DEBUG) {
+          cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
+        }
+      }
+    } else {
+      if (DEBUG) {
+        cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
+      }
+    }
+  } else {
+    if (DEBUG) {
+      cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
+    }
+  }
+
+  if (RecursionLevel < 5) { // Seven is the recursion level threshold.
+    for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+      Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+      if (Walker == Parent) { // don't go back the same bond
+        continue;
+      } else {
+        Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
+      }
+    }
+  }
+  cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
 }
 ;
@@ -1767 +1767 @@
     atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
 {
-        /* OldNormal is normal vector on the old triangle
-         * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
-         * Chord points from b to a!!!
-         */
-        Vector dif_a; //Vector from a to candidate
-        Vector dif_b; //Vector from b to candidate
-        Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
-        Vector TempNormal, Umkreismittelpunkt, Mittelpunkt, helper;
-
-        double CurrentEpsilon = 0.1;
-        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
-        double BallAngle;
-        atom *Walker; // variable atom point
-
-
-        dif_a.CopyVector(&(a->x));
-        dif_a.SubtractVector(&(Candidate->x));
-        dif_b.CopyVector(&(b->x));
-        dif_b.SubtractVector(&(Candidate->x));
-        DirectionCheckPoint.CopyVector(&dif_a);
-        DirectionCheckPoint.Scale(-1);
-        DirectionCheckPoint.ProjectOntoPlane(Chord);
-
-        SideA = dif_b.Norm();
-        SideB = dif_a.Norm();
-        SideC = Chord->Norm();
-        //Chord->Scale(-1);
-
-        alpha = Chord->Angle(&dif_a);
-        beta = M_PI - Chord->Angle(&dif_b);
-        gamma = dif_a.Angle(&dif_b);
-
-
-        if (DEBUG) {
-                cout << "Atom number" << Candidate->nr << endl;
-                Candidate->x.Output((ofstream *) &cout);
-                cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
-        }
-
-        if (a != Candidate and b != Candidate) {
-                //      alpha = dif_a.Angle(&dif_b) / 2.;
-                //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
-                //      SideB = dif_a.Norm();
-                //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
-                //          alpha); // note this is squared of center line length
-                //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-                // Those are remains from Freddie. Needed?
-
-                Umkreisradius = SideA / 2.0 / sin(alpha);
-                //cout << Umkreisradius << endl;
-                //cout << SideB / 2.0 / sin(beta) << endl;
-                //cout << SideC / 2.0 / sin(gamma) << endl;
-
-                if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
-                        // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
-                        Umkreismittelpunkt.Zero();
-                        helper.CopyVector(&a->x);
-                        helper.Scale(sin(2.*alpha));
-                        Umkreismittelpunkt.AddVector(&helper);
-                        helper.CopyVector(&b->x);
-                        helper.Scale(sin(2.*beta));
-                        Umkreismittelpunkt.AddVector(&helper);
-                        helper.CopyVector(&Candidate->x);
-                        helper.Scale(sin(2.*gamma));
-                        Umkreismittelpunkt.AddVector(&helper);
-                        //Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
-                        Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
-                        TempNormal.CopyVector(&dif_a);
-                        TempNormal.VectorProduct(&dif_b);
-                        if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
-                                TempNormal.Scale(-1);
-                        }
-                        TempNormal.Normalize();
-                        Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-                        TempNormal.Scale(Restradius);
-
-                        Mittelpunkt.CopyVector(&Umkreismittelpunkt);
-                        Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
-
-                        AngleCheck.CopyVector(Chord);
-                        AngleCheck.Scale(-0.5);
-                        AngleCheck.SubtractVector(&(b->x));
-                        AngleCheckReference.CopyVector(&AngleCheck);
-                        AngleCheckReference.AddVector(Opt_Mittelpunkt);
-                        AngleCheck.AddVector(&Mittelpunkt);
-
-                        BallAngle = AngleCheck.Angle(&AngleCheckReference);
-
-                        d1->ProjectOntoPlane(&AngleCheckReference);
-                        sign = AngleCheck.ScalarProduct(d1);
-                        if (fabs(sign) < MYEPSILON)
-                                sign = 0;
-                        else
-                                sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-
-
-                        if (Storage[0]< -1.5) { // first Candidate at all
-                                cout << "Next better candidate is " << *Candidate << " with ";
-                                Opt_Candidate = Candidate;
-                                Storage[0] = sign;
-                                Storage[1] = BallAngle;
-                                Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                                cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-                        } else {
-                                /*
-                                 * removed due to change in criterium, now checking angle of ball to old normal.
-                                //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
-                                //within the ball.
-
-                                Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
-                                //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
-
-
-                                if (Distance >RADIUS) { // We have no interference and may now check whether the new point is better.
-                                 */
-                                        //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
-
-                                if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))) { //This will give absolute preference to those in "right-hand" quadrants
-                                        //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
-                                        cout << "Next better candidate is " << *Candidate << " with ";
-                                        Opt_Candidate = Candidate;
-                                        Storage[0] = sign;
-                                        Storage[1] = BallAngle;
-                                        Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                                        cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-                                } else {
-                                        if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0 && Storage[1] > BallAngle) || (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0 && Storage[1] < BallAngle)) {
-                                                //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-                                                cout << "Next better candidate is " << *Candidate << " with ";
-                                                Opt_Candidate = Candidate;
-                                                Storage[0] = sign;
-                                                Storage[1] = BallAngle;
-                                                Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                                                cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-                                        }
-                                }
-                        }
+  /* OldNormal is normal vector on the old triangle
+   * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
+  Vector TempNormal, Umkreismittelpunkt, Mittelpunkt, helper;
+
+  double CurrentEpsilon = 0.1;
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
+  double BallAngle;
+  atom *Walker; // variable atom point
+
+
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
+  DirectionCheckPoint.CopyVector(&dif_a);
+  DirectionCheckPoint.Scale(-1);
+  DirectionCheckPoint.ProjectOntoPlane(Chord);
+
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+
+
+  if (DEBUG) {
+    cout << "Atom number" << Candidate->nr << endl;
+    Candidate->x.Output((ofstream *) &cout);
+    cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
+  }
+
+  if (a != Candidate and b != Candidate) {
+    //      alpha = dif_a.Angle(&dif_b) / 2.;
+    //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
+    //      SideB = dif_a.Norm();
+    //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
+    //          alpha); // note this is squared of center line length
+    //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
+    // Those are remains from Freddie. Needed?
+
+    Umkreisradius = SideA / 2.0 / sin(alpha);
+    //cout << Umkreisradius << endl;
+    //cout << SideB / 2.0 / sin(beta) << endl;
+    //cout << SideC / 2.0 / sin(gamma) << endl;
+
+    if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
+      // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
+      Umkreismittelpunkt.Zero();
+      helper.CopyVector(&a->x);
+      helper.Scale(sin(2.*alpha));
+      Umkreismittelpunkt.AddVector(&helper);
+      helper.CopyVector(&b->x);
+      helper.Scale(sin(2.*beta));
+      Umkreismittelpunkt.AddVector(&helper);
+      helper.CopyVector(&Candidate->x);
+      helper.Scale(sin(2.*gamma));
+      Umkreismittelpunkt.AddVector(&helper);
+      //Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
+      Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+
+      TempNormal.CopyVector(&dif_a);
+      TempNormal.VectorProduct(&dif_b);
+      if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
+        TempNormal.Scale(-1);
+      }
+      TempNormal.Normalize();
+      Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+      TempNormal.Scale(Restradius);
+
+      Mittelpunkt.CopyVector(&Umkreismittelpunkt);
+      Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
+
+      AngleCheck.CopyVector(Chord);
+      AngleCheck.Scale(-0.5);
+      AngleCheck.SubtractVector(&(b->x));
+      AngleCheckReference.CopyVector(&AngleCheck);
+      AngleCheckReference.AddVector(Opt_Mittelpunkt);
+      AngleCheck.AddVector(&Mittelpunkt);
+
+      BallAngle = AngleCheck.Angle(&AngleCheckReference);
+
+      d1->ProjectOntoPlane(&AngleCheckReference);
+      sign = AngleCheck.ScalarProduct(d1);
+      if (fabs(sign) < MYEPSILON)
+        sign = 0;
+      else
+        sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+
+
+      if (Storage[0]< -1.5) { // first Candidate at all
+        cout << "Next better candidate is " << *Candidate << " with ";
+        Opt_Candidate = Candidate;
+        Storage[0] = sign;
+        Storage[1] = BallAngle;
+        Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+        cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+      } else {
+        /*
+         * removed due to change in criterium, now checking angle of ball to old normal.
+        //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
+        //within the ball.
+
+        Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
+        //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+
+
+        if (Distance >RADIUS) { // We have no interference and may now check whether the new point is better.
+         */
+          //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
+
+        if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))) { //This will give absolute preference to those in "right-hand" quadrants
+          //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
+          cout << "Next better candidate is " << *Candidate << " with ";
+          Opt_Candidate = Candidate;
+          Storage[0] = sign;
+          Storage[1] = BallAngle;
+          Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+          cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+        } else {
+          if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0 && Storage[1] > BallAngle) || (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0 && Storage[1] < BallAngle)) {
+            //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+            cout << "Next better candidate is " << *Candidate << " with ";
+            Opt_Candidate = Candidate;
+            Storage[0] = sign;
+            Storage[1] = BallAngle;
+            Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+            cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+          }
+        }
+      }
 /*
                * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
@@ -1947 +1947 @@
                 }
                 */
-                } else {
-                        if (DEBUG) {
-                                cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
-                        }
-                }
-        } else {
-                if (DEBUG) {
-                        cout << "identisch mit Ursprungslinie" << endl;
-                }
-        }
-
-        if (RecursionLevel < 9) { // Five is the recursion level threshold.
-                for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
-                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                        if (Walker == Parent) { // don't go back the same bond
-                                continue;
-                        } else {
-                                Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel+1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS, mol); //call function again
-                        }
-                }
-        }
+    } else {
+      if (DEBUG) {
+        cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+      }
+    }
+  } else {
+    if (DEBUG) {
+      cout << "identisch mit Ursprungslinie" << endl;
+    }
+  }
+
+  if (RecursionLevel < 9) { // Five is the recursion level threshold.
+    for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+      Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+      if (Walker == Parent) { // don't go back the same bond
+        continue;
+      } else {
+        Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel+1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS, mol); //call function again
+      }
+    }
+  }
 };
 
@@ -1984 +1984 @@
     const double& RADIUS, int N, const char *tempbasename)
 {
-        cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
-        Vector direction1;
-        Vector helper;
-        Vector Chord;
-        ofstream *tempstream = NULL;
-        char NumberName[255];
-        double tmp;
-        //atom* Walker;
-        atom* OldThirdPoint;
-
-        double Storage[3];
-        Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
-        Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
-        Storage[2] = 9999999.;
-        atom* Opt_Candidate = NULL;
-        Vector Opt_Mittelpunkt;
-
-        cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-
-        helper.CopyVector(&(Line.endpoints[0]->node->x));
-        for (int i = 0; i < 3; i++) {
-                if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr) {
-                        OldThirdPoint = T.endpoints[i]->node;
-                        helper.SubtractVector(&T.endpoints[i]->node->x);
-                        break;
-                }
-        }
-
-        direction1.CopyVector(&Line.endpoints[0]->node->x);
-        direction1.SubtractVector(&Line.endpoints[1]->node->x);
-        direction1.VectorProduct(&(T.NormalVector));
-
-        if (direction1.ScalarProduct(&helper) < 0) {
-                direction1.Scale(-1);
-        }
-        cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
-
-        Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
-        Chord.SubtractVector(&(Line.endpoints[1]->node->x));
-        cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
-
-
-        Vector Umkreismittelpunkt, a, b, c;
-        double alpha, beta, gamma;
-        a.CopyVector(&(T.endpoints[0]->node->x));
-        b.CopyVector(&(T.endpoints[1]->node->x));
-        c.CopyVector(&(T.endpoints[2]->node->x));
-        a.SubtractVector(&(T.endpoints[1]->node->x));
-        b.SubtractVector(&(T.endpoints[2]->node->x));
-        c.SubtractVector(&(T.endpoints[0]->node->x));
-
-        alpha = M_PI - a.Angle(&c);
-        beta = M_PI - b.Angle(&a);
-        gamma = M_PI - c.Angle(&b);
-        if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-                cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-
-        Umkreismittelpunkt.Zero();
-        helper.CopyVector(&T.endpoints[0]->node->x);
-        helper.Scale(sin(2.*alpha));
-        Umkreismittelpunkt.AddVector(&helper);
-        helper.CopyVector(&T.endpoints[1]->node->x);
-        helper.Scale(sin(2.*beta));
-        Umkreismittelpunkt.AddVector(&helper);
-        helper.CopyVector(&T.endpoints[2]->node->x);
-        helper.Scale(sin(2.*gamma));
-        Umkreismittelpunkt.AddVector(&helper);
-        //Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
-        //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-        Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-        //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-        cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
-        cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
-
-        cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-        cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
-        if (DEBUG)
-                cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
-        tmp = 0;
-        for (int i=0;i<NDIM;i++) {
-                helper.CopyVector(&T.endpoints[i]->node->x);
-                helper.SubtractVector(&Umkreismittelpunkt);
-                if (DEBUG)
-                        cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
-                if (tmp == 0) // set on first time for comparison against next ones
-                        tmp = helper.Norm();
-                if (fabs(helper.Norm() - tmp) > MYEPSILON)
-                        cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
-        }
-
-        cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
-
-        Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-        Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-        if (Opt_Candidate == NULL) {
-                cerr << "WARNING: Could not find a suitable candidate." << endl;
-                return false;
-        }
-        cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << endl;
-
-        // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-        bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node);
-
-        if (flag) { // if so, add
-                AddTrianglePoint(Opt_Candidate, 0);
-                AddTrianglePoint(Line.endpoints[0]->node, 1);
-                AddTrianglePoint(Line.endpoints[1]->node, 2);
-
-                AddTriangleLine(TPS[0], TPS[1], 0);
-                AddTriangleLine(TPS[0], TPS[2], 1);
-                AddTriangleLine(TPS[1], TPS[2], 2);
-
-                BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-                AddTriangleToLines();
-
-                BTS->GetNormalVector(BTS->NormalVector);
-
-                if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  || (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) || (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) ) {
-                        BTS->NormalVector.Scale(-1);
-                };
-                cout << Verbose(1) << "New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-                cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
-        } else { // else, yell and do nothing
-                cout << Verbose(1) << "This triangle consisting of ";
-                cout << *Opt_Candidate << ", ";
-                cout << *Line.endpoints[0]->node << " and ";
-                cout << *Line.endpoints[1]->node << " ";
-                cout << "is invalid!" << endl;
-                return false;
-        }
-
-        if ((TrianglesOnBoundaryCount % 10) == 0) {
-                sprintf(NumberName, "-%d", TriangleFilesWritten);
-                if (DoTecplotOutput) {
-                        string NameofTempFile(tempbasename);
-                        NameofTempFile.append(NumberName);
-                        NameofTempFile.append(TecplotSuffix);
-                        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-                        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-                        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
-                        tempstream->close();
-                        tempstream->flush();
-                        delete(tempstream);
-                }
-
-                if (DoRaster3DOutput) {
-                        string NameofTempFile(tempbasename);
-                        NameofTempFile.append(NumberName);
-                        NameofTempFile.append(Raster3DSuffix);
-                        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-                        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-                        write_raster3d_file(out, tempstream, this, mol);
-                        tempstream->close();
-                        tempstream->flush();
-                        delete(tempstream);
-                }
-                if (DoTecplotOutput || DoRaster3DOutput)
-                        TriangleFilesWritten++;
-        }
-
-        cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
-        return true;
+  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+  Vector direction1;
+  Vector helper;
+  Vector Chord;
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  double tmp;
+  //atom* Walker;
+  atom* OldThirdPoint;
+
+  double Storage[3];
+  Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[2] = 9999999.;
+  atom* Opt_Candidate = NULL;
+  Vector Opt_Mittelpunkt;
+
+  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+
+  helper.CopyVector(&(Line.endpoints[0]->node->x));
+  for (int i = 0; i < 3; i++) {
+    if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr) {
+      OldThirdPoint = T.endpoints[i]->node;
+      helper.SubtractVector(&T.endpoints[i]->node->x);
+      break;
+    }
+  }
+
+  direction1.CopyVector(&Line.endpoints[0]->node->x);
+  direction1.SubtractVector(&Line.endpoints[1]->node->x);
+  direction1.VectorProduct(&(T.NormalVector));
+
+  if (direction1.ScalarProduct(&helper) < 0) {
+    direction1.Scale(-1);
+  }
+  cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
+
+  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
+  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
+  cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
+
+
+  Vector Umkreismittelpunkt, a, b, c;
+  double alpha, beta, gamma;
+  a.CopyVector(&(T.endpoints[0]->node->x));
+  b.CopyVector(&(T.endpoints[1]->node->x));
+  c.CopyVector(&(T.endpoints[2]->node->x));
+  a.SubtractVector(&(T.endpoints[1]->node->x));
+  b.SubtractVector(&(T.endpoints[2]->node->x));
+  c.SubtractVector(&(T.endpoints[0]->node->x));
+
+  alpha = M_PI - a.Angle(&c);
+  beta = M_PI - b.Angle(&a);
+  gamma = M_PI - c.Angle(&b);
+  if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
+    cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+
+  Umkreismittelpunkt.Zero();
+  helper.CopyVector(&T.endpoints[0]->node->x);
+  helper.Scale(sin(2.*alpha));
+  Umkreismittelpunkt.AddVector(&helper);
+  helper.CopyVector(&T.endpoints[1]->node->x);
+  helper.Scale(sin(2.*beta));
+  Umkreismittelpunkt.AddVector(&helper);
+  helper.CopyVector(&T.endpoints[2]->node->x);
+  helper.Scale(sin(2.*gamma));
+  Umkreismittelpunkt.AddVector(&helper);
+  //Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+  //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+  Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+  //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+  cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
+  cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
+
+  cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+  cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
+  if (DEBUG)
+    cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
+  tmp = 0;
+  for (int i=0;i<NDIM;i++) {
+    helper.CopyVector(&T.endpoints[i]->node->x);
+    helper.SubtractVector(&Umkreismittelpunkt);
+    if (DEBUG)
+      cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
+    if (tmp == 0) // set on first time for comparison against next ones
+      tmp = helper.Norm();
+    if (fabs(helper.Norm() - tmp) > MYEPSILON)
+      cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
+  }
+
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+  if (Opt_Candidate == NULL) {
+    cerr << "WARNING: Could not find a suitable candidate." << endl;
+    return false;
+  }
+  cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << endl;
+
+  // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+  bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node);
+
+  if (flag) { // if so, add
+    AddTrianglePoint(Opt_Candidate, 0);
+    AddTrianglePoint(Line.endpoints[0]->node, 1);
+    AddTrianglePoint(Line.endpoints[1]->node, 2);
+
+    AddTriangleLine(TPS[0], TPS[1], 0);
+    AddTriangleLine(TPS[0], TPS[2], 1);
+    AddTriangleLine(TPS[1], TPS[2], 2);
+
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    AddTriangleToLines();
+
+    BTS->GetNormalVector(BTS->NormalVector);
+
+    if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  || (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) || (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) ) {
+      BTS->NormalVector.Scale(-1);
+    };
+    cout << Verbose(1) << "New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
+    cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
+  } else { // else, yell and do nothing
+    cout << Verbose(1) << "This triangle consisting of ";
+    cout << *Opt_Candidate << ", ";
+    cout << *Line.endpoints[0]->node << " and ";
+    cout << *Line.endpoints[1]->node << " ";
+    cout << "is invalid!" << endl;
+    return false;
+  }
+
+  if ((TrianglesOnBoundaryCount % 10) == 0) {
+    sprintf(NumberName, "-%d", TriangleFilesWritten);
+    if (DoTecplotOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(TecplotSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+
+    if (DoRaster3DOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(Raster3DSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_raster3d_file(out, tempstream, this, mol);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoTecplotOutput || DoRaster3DOutput)
+      TriangleFilesWritten++;
+  }
+
+  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+  return true;
 };
 
@@ -2159 +2159 @@
  */
 bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c) {
-        LineMap::iterator TryAndError;
-        PointTestPair InsertPoint;
-        bool Present[3];
-
-        *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
-        TPS[0] = new class BoundaryPointSet(a);
-        TPS[1] = new class BoundaryPointSet(b);
-        TPS[2] = new class BoundaryPointSet(c);
-        for (int i=0;i<3;i++) { // check through all endpoints
-                InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i]));
-                Present[i] = !InsertPoint.second;
-                if (Present[i]) { // if new point was not present before, increase counter
-                        delete TPS[i];
-                        *out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl;
-                        TPS[i] = (InsertPoint.first)->second;
-                }
-        }
-
-        // check lines
-        for (int i=0;i<3;i++)
-                if (Present[i])
-                        for (int j=i;j<3;j++)
-                                if (Present[j]) {
-                                        TryAndError = TPS[i]->lines.find(TPS[j]->node->nr);
-                                        if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) {
-                                                *out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl;
-                                                *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
-                                                return false;
-                                        }
-                                }
-        *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
-        return true;
+  LineMap::iterator TryAndError;
+  PointTestPair InsertPoint;
+  bool Present[3];
+
+  *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+  TPS[0] = new class BoundaryPointSet(a);
+  TPS[1] = new class BoundaryPointSet(b);
+  TPS[2] = new class BoundaryPointSet(c);
+  for (int i=0;i<3;i++) { // check through all endpoints
+    InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i]));
+    Present[i] = !InsertPoint.second;
+    if (Present[i]) { // if new point was not present before, increase counter
+      delete TPS[i];
+      *out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl;
+      TPS[i] = (InsertPoint.first)->second;
+    }
+  }
+
+  // check lines
+  for (int i=0;i<3;i++)
+    if (Present[i])
+      for (int j=i;j<3;j++)
+        if (Present[j]) {
+          TryAndError = TPS[i]->lines.find(TPS[j]->node->nr);
+          if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) {
+            *out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl;
+            *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+            return false;
+          }
+        }
+  *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+  return true;
 };
 
@@ -2198 +2198 @@
     molecule* mol, double RADIUS)
 {
-        cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
-        int i;
-        Vector AngleCheck;
-        atom* Walker;
-        double norm = -1., angle;
-
-        // check if we only have one unique point yet ...
-        if (a != Candidate) {
-                cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
-                AngleCheck.CopyVector(&(Candidate->x));
-                AngleCheck.SubtractVector(&(a->x));
-                norm = AngleCheck.Norm();
-                // second point shall have smallest angle with respect to Oben vector
-                if (norm < RADIUS) {
-                        angle = AngleCheck.Angle(&Oben);
-                        if (angle < Storage[0]) {
-                                //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                                cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
-                                Opt_Candidate = Candidate;
-                                Storage[0] = AngleCheck.Angle(&Oben);
-                                //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-                        } else {
-                                cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
-                        }
-                } else {
-                        cout << "Refused due to Radius " << norm << endl;
-                }
-        }
-
-        // if not recursed to deeply, look at all its bonds
-        if (RecursionLevel < 7) {
-                for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) {
-                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                        if (Walker == Parent) // don't go back along the bond we came from
-                                continue;
-                        else
-                                Find_second_point_for_Tesselation(a, Walker, Candidate, RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
-                }
-        }
-        cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
+  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
+  int i;
+  Vector AngleCheck;
+  atom* Walker;
+  double norm = -1., angle;
+
+  // check if we only have one unique point yet ...
+  if (a != Candidate) {
+    cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
+    AngleCheck.CopyVector(&(Candidate->x));
+    AngleCheck.SubtractVector(&(a->x));
+    norm = AngleCheck.Norm();
+    // second point shall have smallest angle with respect to Oben vector
+    if (norm < RADIUS) {
+      angle = AngleCheck.Angle(&Oben);
+      if (angle < Storage[0]) {
+        //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+        cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
+        Opt_Candidate = Candidate;
+        Storage[0] = AngleCheck.Angle(&Oben);
+        //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+      } else {
+        cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+      }
+    } else {
+      cout << "Refused due to Radius " << norm << endl;
+    }
+  }
+
+  // if not recursed to deeply, look at all its bonds
+  if (RecursionLevel < 7) {
+    for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) {
+      Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+      if (Walker == Parent) // don't go back along the bond we came from
+        continue;
+      else
+        Find_second_point_for_Tesselation(a, Walker, Candidate, RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+    }
+  }
+  cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
 };
 
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-        cout << Verbose(1) << "Begin of Find_starting_triangle\n";
-        int i = 0;
-        atom* Walker;
-        atom* FirstPoint;
-        atom* SecondPoint;
-        atom* max_index[NDIM];
-        double max_coordinate[NDIM];
-        Vector Oben;
-        Vector helper;
-        Vector Chord;
-        Vector CenterOfFirstLine;
-
-        Oben.Zero();
-
-        for (i = 0; i < 3; i++) {
-                max_index[i] = NULL;
-                max_coordinate[i] = -1;
-        }
-        cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
-
-        // 1. searching topmost atom with respect to each axis
-        Walker = mol->start;
-        while (Walker->next != mol->end) {
-                Walker = Walker->next;
-                for (i = 0; i < 3; i++) {
-                        if (Walker->x.x[i] > max_coordinate[i]) {
-                                max_coordinate[i] = Walker->x.x[i];
-                                max_index[i] = Walker;
-                        }
-                }
-        }
-
-        cout << Verbose(2) << "Found maximum coordinates: ";
-        for (int i=0;i<NDIM;i++)
-                cout << i << ": " << *max_index[i] << "\t";
-        cout << endl;
+  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+  int i = 0;
+  atom* Walker;
+  atom* FirstPoint;
+  atom* SecondPoint;
+  atom* max_index[NDIM];
+  double max_coordinate[NDIM];
+  Vector Oben;
+  Vector helper;
+  Vector Chord;
+  Vector CenterOfFirstLine;
+
+  Oben.Zero();
+
+  for (i = 0; i < 3; i++) {
+    max_index[i] = NULL;
+    max_coordinate[i] = -1;
+  }
+  cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
+
+  // 1. searching topmost atom with respect to each axis
+  Walker = mol->start;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    for (i = 0; i < 3; i++) {
+      if (Walker->x.x[i] > max_coordinate[i]) {
+        max_coordinate[i] = Walker->x.x[i];
+        max_index[i] = Walker;
+      }
+    }
+  }
+
+  cout << Verbose(2) << "Found maximum coordinates: ";
+  for (int i=0;i<NDIM;i++)
+    cout << i << ": " << *max_index[i] << "\t";
+  cout << endl;
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-        const int k = 1;
-        Oben.x[k] = 1.;
-        FirstPoint = max_index[k];
-
-        cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
-        double Storage[3];
-        atom* Opt_Candidate = NULL;
-        Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-        Storage[1] = 999999.; // This will be an angle looking for the third point.
-        Storage[2] = 999999.;
-
-        Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
-        SecondPoint = Opt_Candidate;
-        cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-
-        helper.CopyVector(&(FirstPoint->x));
-        helper.SubtractVector(&(SecondPoint->x));
-        helper.Normalize();
-        Oben.ProjectOntoPlane(&helper);
-        Oben.Normalize();
-        helper.VectorProduct(&Oben);
-        Storage[0] = -2.; // This will indicate the quadrant.
-        Storage[1] = 9999999.; // This will be an angle looking for the third point.
-        Storage[2] = 9999999.;
-
-        Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-        Chord.SubtractVector(&(SecondPoint->x));
-        // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-
-        cout << Verbose(2) << "Looking for third point candidates \n";
-        // look in one direction of baseline for initial candidate
-        Opt_Candidate = NULL;
-        CenterOfFirstLine.CopyVector(&Chord);
-        CenterOfFirstLine.Scale(0.5);
-        CenterOfFirstLine.AddVector(&(SecondPoint->x));
-
-        cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
-        Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
-        // look in other direction of baseline for possible better candidate
-        cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
-        Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
-        cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
-
-        // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
-
-        // Finally, we only have to add the found points
-        AddTrianglePoint(FirstPoint, 0);
-        AddTrianglePoint(SecondPoint, 1);
-        AddTrianglePoint(Opt_Candidate, 2);
-        // ... and respective lines
-        AddTriangleLine(TPS[0], TPS[1], 0);
-        AddTriangleLine(TPS[1], TPS[2], 1);
-        AddTriangleLine(TPS[0], TPS[2], 2);
-        // ... and triangles to the Maps of the Tesselation class
-        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-        AddTriangleToLines();
-        // ... and calculate its normal vector (with correct orientation)
-        Oben.Scale(-1.);
-        BTS->GetNormalVector(Oben);
-        cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
-        cout << Verbose(1) << "End of Find_starting_triangle\n";
+  const int k = 1;
+  Oben.x[k] = 1.;
+  FirstPoint = max_index[k];
+
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
+  double Storage[3];
+  atom* Opt_Candidate = NULL;
+  Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[1] = 999999.; // This will be an angle looking for the third point.
+  Storage[2] = 999999.;
+
+  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
+  SecondPoint = Opt_Candidate;
+  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+
+  helper.CopyVector(&(FirstPoint->x));
+  helper.SubtractVector(&(SecondPoint->x));
+  helper.Normalize();
+  Oben.ProjectOntoPlane(&helper);
+  Oben.Normalize();
+  helper.VectorProduct(&Oben);
+  Storage[0] = -2.; // This will indicate the quadrant.
+  Storage[1] = 9999999.; // This will be an angle looking for the third point.
+  Storage[2] = 9999999.;
+
+  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+  Chord.SubtractVector(&(SecondPoint->x));
+  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+
+  cout << Verbose(2) << "Looking for third point candidates \n";
+  // look in one direction of baseline for initial candidate
+  Opt_Candidate = NULL;
+  CenterOfFirstLine.CopyVector(&Chord);
+  CenterOfFirstLine.Scale(0.5);
+  CenterOfFirstLine.AddVector(&(SecondPoint->x));
+
+  cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
+  // look in other direction of baseline for possible better candidate
+  cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
+  cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+
+  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+
+  // Finally, we only have to add the found points
+  AddTrianglePoint(FirstPoint, 0);
+  AddTrianglePoint(SecondPoint, 1);
+  AddTrianglePoint(Opt_Candidate, 2);
+  // ... and respective lines
+  AddTriangleLine(TPS[0], TPS[1], 0);
+  AddTriangleLine(TPS[1], TPS[2], 1);
+  AddTriangleLine(TPS[0], TPS[2], 2);
+  // ... and triangles to the Maps of the Tesselation class
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTriangleToLines();
+  // ... and calculate its normal vector (with correct orientation)
+  Oben.Scale(-1.);
+  BTS->GetNormalVector(Oben);
+  cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
+  cout << Verbose(1) << "End of Find_starting_triangle\n";
 };
 
 void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol, const char *filename, const double RADIUS)
 {
-        int N = 0;
-        struct Tesselation *Tess = new Tesselation;
-        cout << Verbose(1) << "Entering search for non convex hull. " << endl;
-        cout << flush;
-        LineMap::iterator baseline;
-        cout << Verbose(0) << "Begin of Find_non_convex_border\n";
-        bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
-        bool failflag = false;
-        if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
-                mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
-
-        Tess->Find_starting_triangle(mol, RADIUS);
-
-        baseline = Tess->LinesOnBoundary.begin();
-        while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
-                if (baseline->second->TrianglesCount == 1) {
-                        failflag = Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
-                        flag = flag || failflag;
-                        if (!failflag)
-                                cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
-                } else {
-                        cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
-                }
-                N++;
-                baseline++;
-                if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
-                        baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-                        flag = false;
-                }
-        }
-        if (failflag) {
-                cout << Verbose(1) << "Writing final tecplot file\n";
-                if (DoTecplotOutput)
-                        write_tecplot_file(out, tecplot, Tess, mol, -1);
-                if (DoRaster3DOutput)
-                        write_raster3d_file(out, tecplot, Tess, mol);
-        } else {
-                cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
-        }
-
-        cout << Verbose(0) << "End of Find_non_convex_border\n";
-        delete(Tess);
+  int N = 0;
+  struct Tesselation *Tess = new Tesselation;
+  cout << Verbose(1) << "Entering search for non convex hull. " << endl;
+  cout << flush;
+  LineMap::iterator baseline;
+  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
+  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
+  bool failflag = false;
+  if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
+    mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
+
+  Tess->Find_starting_triangle(mol, RADIUS);
+
+  baseline = Tess->LinesOnBoundary.begin();
+  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
+    if (baseline->second->TrianglesCount == 1) {
+      failflag = Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
+      flag = flag || failflag;
+      if (!failflag)
+        cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
+    } else {
+      cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
+    }
+    N++;
+    baseline++;
+    if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
+      baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      flag = false;
+    }
+  }
+  if (failflag) {
+    cout << Verbose(1) << "Writing final tecplot file\n";
+    if (DoTecplotOutput)
+      write_tecplot_file(out, tecplot, Tess, mol, -1);
+    if (DoRaster3DOutput)
+      write_raster3d_file(out, tecplot, Tess, mol);
+  } else {
+    cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
+  }
+
+  cout << Verbose(0) << "End of Find_non_convex_border\n";
+  delete(Tess);
 };
 

src/builder.cpp

@@ -562 +562 @@
       break;
     case 'd':
-        cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
-        cout << Verbose(0) << "Shall we rotate? [0/1]: ";
-        cin >> Z;
-        if ((Z >=0) && (Z <=1))
-          mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
-        else
-          mol->PrincipalAxisSystem((ofstream *)&cout, false);
-        break;
+      cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+      cout << Verbose(0) << "Shall we rotate? [0/1]: ";
+      cin >> Z;
+      if ((Z >=0) && (Z <=1))
+        mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
+      else
+        mol->PrincipalAxisSystem((ofstream *)&cout, false);
+      break;
     case 'e':
-        cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-        VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
-        break;
+      cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+      VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
+      break;
     case 'f':
       mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
@@ -1020 +1020 @@
               }
               break;
-                                                case 'A':
-                                                        ExitFlag = 1;
-                                                        if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-')) {
-                                                                ExitFlag =255;
-                                                                cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
-                                                        } else {
-                                                                cout << "Parsing bonds from " << argv[argptr] << "." << endl;
-                                                                ifstream *input = new ifstream(argv[argptr]);
-                                                                mol->CreateAdjacencyList2((ofstream *)&cout, input);
-                                                                input->close();
-                                                        }
-                                                        break;
-                                                case 'N':
-                                                        ExitFlag = 1;
-                                                        if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
-                                                                ExitFlag = 255;
-                                                                cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
-                                                        } else {
-                                                                cout << Verbose(0) << "Evaluating npn-convex envelope.";
-                                                                string TempName(argv[argptr+1]);
-                                                                TempName.append(".r3d");
-                                                                ofstream *output = new ofstream(TempName.c_str(), ios::trunc);
-                                                                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
-                                                                Find_non_convex_border((ofstream *)&cout, output, mol, argv[argptr+1], atof(argv[argptr]));
-                                                                output->close();
-                                                                delete(output);
-                                                                argptr+=2;
-                                                        }
-                                                        break;
+            case 'A':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-')) {
+                ExitFlag =255;
+                cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
+              } else {
+                cout << "Parsing bonds from " << argv[argptr] << "." << endl;
+                ifstream *input = new ifstream(argv[argptr]);
+                mol->CreateAdjacencyList2((ofstream *)&cout, input);
+                input->close();
+              }
+              break;
+            case 'N':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
+              } else {
+                cout << Verbose(0) << "Evaluating npn-convex envelope.";
+                string TempName(argv[argptr+1]);
+                TempName.append(".r3d");
+                ofstream *output = new ofstream(TempName.c_str(), ios::trunc);
+                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
+                Find_non_convex_border((ofstream *)&cout, output, mol, argv[argptr+1], atof(argv[argptr]));
+                output->close();
+                delete(output);
+                argptr+=2;
+              }
+              break;
             case 'T':
               ExitFlag = 1;
@@ -1369 +1369 @@
   char choice;  // menu choice char
   Vector x,y,z,n;  // coordinates for absolute point in cell volume
-        double *factor; // unit factor if desired
+  double *factor;  // unit factor if desired
   bool valid; // flag if input was valid or not
   ifstream test;

src/config.cpp

@@ -1082 +1082 @@
 {
   bool result = true;
-        // bring MaxTypes up to date
-        mol->CountElements();
+  // bring MaxTypes up to date
+  mol->CountElements();
   ofstream *output = NULL;
   output = new ofstream(filename, ios::out);

src/datacreator.cpp

@@ -335 +335 @@
     for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
       for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
-        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
-        norm += tmp*tmp;
+        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
+        norm += tmp*tmp;
       }
     }

src/helpers.cpp

@@ -73 +73 @@
   void *buffer = NULL;
   if (OldPointer == NULL)
-        //cout << Verbose(0) << "ReAlloc impossible - old is NULL: " << output << endl;
+    //cout << Verbose(0) << "ReAlloc impossible - old is NULL: " << output << endl;
     buffer = (void *)malloc(size); // malloc
   else

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--;) {

src/parser.cpp

@@ -786 +786 @@
           }
           if (k != -1) {
-                //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
+            //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
             Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
           }
@@ -855 +855 @@
                   }
                   if (Order == SubOrder) { // equal order is always copy from Energies
-                        //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
                   } else {
-                        //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
                   }

src/vector.cpp

@@ -327 +327 @@
   }
   ost << ")";
-        return ost;
+  return ost;
 };
 
@@ -517 +517 @@
 bool Vector::MakeNormalVector(const Vector *y1)
 {
-        bool result = false;
+  bool result = false;
   Vector x1;
   x1.CopyVector(y1);
@@ -523 +523 @@
   SubtractVector(&x1);
   for (int i=NDIM;i--;)
-          result = result || (fabs(x[i]) > MYEPSILON);
+    result = result || (fabs(x[i]) > MYEPSILON);
 
   return result;
@@ -597 +597 @@
 bool Vector::LSQdistance(Vector **vectors, int num)
 {
-        int j;
-
-        for (j=0;j<num;j++) {
-                cout << Verbose(1) << j << "th atom's vector: ";
-                (vectors[j])->Output((ofstream *)&cout);
-                cout << endl;
-        }
+  int j;
+
+  for (j=0;j<num;j++) {
+    cout << Verbose(1) << j << "th atom's vector: ";
+    (vectors[j])->Output((ofstream *)&cout);
+    cout << endl;
+  }
 
   int np = 3;
-        struct LSQ_params par;
+  struct LSQ_params par;
 
    const gsl_multimin_fminimizer_type *T =
@@ -627 +627 @@
    /* Starting point */
    par.vectors = vectors;
-         par.num = num;
-
-         for (i=NDIM;i--;)
-                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
+   par.num = num;
+
+   for (i=NDIM;i--;)
+    gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
 
    /* Initialize method and iterate */
@@ -671 +671 @@
    gsl_multimin_fminimizer_free (s);
 
-        return true;
+  return true;
 };
 

src/verbose.cpp

@@ -22 +22 @@
 ostream& operator<<(ostream& ost,const Verbose& m)
 {
-        return m.print(ost);
+  return m.print(ost);
 };
 
@@ -33 +33 @@
 ostream& Binary::print (ostream &ost) const
 {
-        int bits = 1, counter = 1;
-        while ((bits = 1 << counter) < BinaryNumber)
+  int bits = 1, counter = 1;
+  while ((bits = 1 << counter) < BinaryNumber)
     counter++;
   for (int i=0;i<counter-1;i++) {
-        if ((BinaryNumber & (1 << i)) == 0)
-                ost.put('0');
-                else
-                ost.put('1');
+    if ((BinaryNumber & (1 << i)) == 0)
+      ost.put('0');
+     else
+      ost.put('1');
   }
   ost.put('b');
@@ -56 +56 @@
 ostream& operator<<(ostream& ost,const Binary& m)
 {
-        return m.print(ost);
+  return m.print(ost);
 };