Changeset cee0b57 for src


Ignore:
Timestamp:
Oct 5, 2009, 10:10:53 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:
ebcade
Parents:
d09ff7
Message:

class molecule implementation split up into six separate parts.

  • dynamics: Verlet integration and constraint potential
  • fragmentation: BOSSANOVA scheme
  • geometry: all operations acting on the Vector's inside the atom's
  • graph: supplementary functions for fragmentation, treating molecule as a bonding graph
  • pointcloud: implementations of virtual functions for pointcloud class, needed for Tesselation
Location:
src
Files:
5 added
10 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    rd09ff7 rcee0b57  
    1 SOURCE = atom.cpp bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
    2 HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp memoryallocator.hpp memoryusageobserver.cpp memoryusageobserver.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
     1SOURCE = atom.cpp bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
     2HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp memoryallocator.hpp memoryusageobserver.cpp memoryusageobserver.hpp molecule.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
    33
    44noinst_PROGRAMS =  ActOnAllTest VectorUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest TesselationUnitTest

  • src/boundary.hpp

    rd09ff7 rcee0b57  
    1212#include "config.hpp"
    1313#include "linkedcell.hpp"
    14 #include "molecules.hpp"
     14#include "molecule.hpp"
    1515#include "tesselation.hpp"
    1616

  • src/builder.cpp

    rd09ff7 rcee0b57  
    5454#include "helpers.hpp"
    5555#include "memoryusageobserverunittest.hpp"
    56 #include "molecules.hpp"
     56#include "molecule.hpp"
    5757/********************************************* Subsubmenu routine ************************************/
    5858

  • src/config.hpp

    rd09ff7 rcee0b57  
    1616#endif
    1717
    18 #include "molecules.hpp"
     18#include "molecule.hpp"
    1919#include "periodentafel.hpp"
    2020

  • src/graph.hpp

    rd09ff7 rcee0b57  
    1717#include <multimap>
    1818
    19 #include "molecules.hpp"
     19#include "molecule.hpp"
    2020
    2121class Graph;

  • src/linkedcell.cpp

    rd09ff7 rcee0b57  
    77
    88#include "linkedcell.hpp"
    9 #include "molecules.hpp"
     9#include "molecule.hpp"
    1010#include "tesselation.hpp"
    1111

  • src/molecule.cpp

    rd09ff7 rcee0b57  
    77#include "config.hpp"
    88#include "memoryallocator.hpp"
    9 #include "molecules.hpp"
     9#include "molecule.hpp"
    1010
    1111/************************************* Functions for class molecule *********************************/
    … …  
    6868};
    6969
    70 
    71 /** Determine center of all atoms.
    72  * \param *out output stream for debugging
    73  * \return pointer to allocated with central coordinates
    74  */
    75 Vector *molecule::GetCenter(ofstream *out)
    76 {
    77   Vector *center = DetermineCenterOfAll(out);
    78   return center;
    79 };
    80 
    81 /** Return current atom in the list.
    82  * \return pointer to atom or NULL if none present
    83  */
    84 TesselPoint *molecule::GetPoint()
    85 {
    86   if ((InternalPointer != start) && (InternalPointer != end))
    87     return InternalPointer;
    88   else
    89     return NULL;
    90 };
    91 
    92 /** Return pointer to one after last atom in the list.
    93  * \return pointer to end marker
    94  */
    95 TesselPoint *molecule::GetTerminalPoint()
    96 {
    97   return end;
    98 };
    99 
    100 /** Go to next atom.
    101  * Stops at last one.
    102  */
    103 void molecule::GoToNext()
    104 {
    105   if (InternalPointer != end)
    106     InternalPointer = InternalPointer->next;
    107 };
    108 
    109 /** Go to previous atom.
    110  * Stops at first one.
    111  */
    112 void molecule::GoToPrevious()
    113 {
    114   if (InternalPointer->previous != start)
    115     InternalPointer = InternalPointer->previous;
    116 };
    117 
    118 /** Goes to first atom.
    119  */
    120 void molecule::GoToFirst()
    121 {
    122   InternalPointer = start->next;
    123 };
    124 
    125 /** Goes to last atom.
    126  */
    127 void molecule::GoToLast()
    128 {
    129   InternalPointer = end->previous;
    130 };
    131 
    132 /** Checks whether we have any atoms in molecule.
    133  * \return true - no atoms, false - not empty
    134  */
    135 bool molecule::IsEmpty()
    136 {
    137   return (start->next == end);
    138 };
    139 
    140 /** Checks whether we are at the last atom
    141  * \return true - current atom is last one, false - is not last one
    142  */
    143 bool molecule::IsEnd()
    144 {
    145   return (InternalPointer == end);
    146 };
    14770
    14871/** Adds given atom \a *pointer from molecule list.
    … …  
    575498  while(Binder->next != last) {
    576499    Binder = Binder->next;
     500
    577501    // get the pendant atoms of current bond in the copy molecule
    578     LeftAtom = copy->start;
    579     while (LeftAtom->next != copy->end) {
    580       LeftAtom = LeftAtom->next;
    581       if (LeftAtom->father == Binder->leftatom)
    582         break;
    583     }
    584     RightAtom = copy->start;
    585     while (RightAtom->next != copy->end) {
    586       RightAtom = RightAtom->next;
    587       if (RightAtom->father == Binder->rightatom)
    588         break;
    589     }
     502    copy->ActOnAllAtoms( &atom::EqualsFather, Binder->leftatom, &LeftAtom );
     503    copy->ActOnAllAtoms( &atom::EqualsFather, Binder->rightatom, &RightAtom );
     504
    590505    NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
    591506    NewBond->Cyclic = Binder->Cyclic;
    … …  
    595510  }
    596511  // correct fathers
    597   Walker = copy->start;
    598   while(Walker->next != copy->end) {
    599     Walker = Walker->next;
    600     if (Walker->father->father == Walker->father)   // same atom in copy's father points to itself
    601       Walker->father = Walker;  // set father to itself (copy of a whole molecule)
    602     else
    603      Walker->father = Walker->father->father;  // set father to original's father
    604   }
     512  ActOnAllAtoms( &atom::CorrectFather );
     513
    605514  // copy values
    606515  copy->CountAtoms((ofstream *)&cout);
    … …  
    644553 * \return pointer to bond or NULL on failure
    645554 */
    646 bond * molecule::AddBond(atom *atom1, atom *atom2, int degree=1)
     555bond * molecule::AddBond(atom *atom1, atom *atom2, int degree)
    647556{
    648557  bond *Binder = NULL;
    … …  
    714623};
    715624
    716 /** Centers the molecule in the box whose lengths are defined by vector \a *BoxLengths.
    717  * \param *out output stream for debugging
    718  */
    719 bool molecule::CenterInBox(ofstream *out)
    720 {
    721   bool status = true;
    722   Vector x;
    723   double *M = ReturnFullMatrixforSymmetric(cell_size);
    724   double *Minv = x.InverseMatrix(M);
    725   Vector *Center = DetermineCenterOfAll(out);
    726 
    727   // go through all atoms
    728   atom *ptr = start;  // start at first in list
    729   while (ptr->next != end) {  // continue with second if present
    730     ptr = ptr->next;
    731     //ptr->Output(1,1,out);
    732     // multiply its vector with matrix inverse
    733     x.CopyVector(&ptr->x);
    734     x.SubtractVector(Center); // now, it's centered at origin
    735     x.MatrixMultiplication(Minv);
    736     // truncate to [0,1] for each axis
    737     for (int i=0;i<NDIM;i++) {
    738       x.x[i] += 0.5;  // set to center of box
    739       while (x.x[i] >= 1.)
    740         x.x[i] -= 1.;
    741       while (x.x[i] < 0.)
    742         x.x[i] += 1.;
    743     }
    744     x.MatrixMultiplication(M);
    745     ptr->x.CopyVector(&x);
    746   }
    747   delete(M);
    748   delete(Minv);
    749   delete(Center);
    750   return status;
    751 };
    752 
    753 
    754 /** Bounds the molecule in the box whose lengths are defined by vector \a *BoxLengths.
    755  * \param *out output stream for debugging
    756  */
    757 bool molecule::BoundInBox(ofstream *out)
    758 {
    759   bool status = true;
    760   Vector x;
    761   double *M = ReturnFullMatrixforSymmetric(cell_size);
    762   double *Minv = x.InverseMatrix(M);
    763 
    764   // go through all atoms
    765   atom *ptr = start;  // start at first in list
    766   while (ptr->next != end) {  // continue with second if present
    767     ptr = ptr->next;
    768     //ptr->Output(1,1,out);
    769     // multiply its vector with matrix inverse
    770     x.CopyVector(&ptr->x);
    771     x.MatrixMultiplication(Minv);
    772     // truncate to [0,1] for each axis
    773     for (int i=0;i<NDIM;i++) {
    774       while (x.x[i] >= 1.)
    775         x.x[i] -= 1.;
    776       while (x.x[i] < 0.)
    777         x.x[i] += 1.;
    778     }
    779     x.MatrixMultiplication(M);
    780     ptr->x.CopyVector(&x);
    781   }
    782   delete(M);
    783   delete(Minv);
    784   return status;
    785 };
    786 
    787 /** Centers the edge of the atoms at (0,0,0).
    788  * \param *out output stream for debugging
    789  * \param *max coordinates of other edge, specifying box dimensions.
    790  */
    791 void molecule::CenterEdge(ofstream *out, Vector *max)
    792 {
    793   Vector *min = new Vector;
    794 
    795 //  *out << Verbose(3) << "Begin of CenterEdge." << endl;
    796   atom *ptr = start->next;  // start at first in list
    797   if (ptr != end) {   //list not empty?
    798     for (int i=NDIM;i--;) {
    799       max->x[i] = ptr->x.x[i];
    800       min->x[i] = ptr->x.x[i];
    801     }
    802     while (ptr->next != end) {  // continue with second if present
    803       ptr = ptr->next;
    804       //ptr->Output(1,1,out);
    805       for (int i=NDIM;i--;) {
    806         max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
    807         min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i];
    808       }
    809     }
    810 //    *out << Verbose(4) << "Maximum is ";
    811 //    max->Output(out);
    812 //    *out << ", Minimum is ";
    813 //    min->Output(out);
    814 //    *out << endl;
    815     min->Scale(-1.);
    816     max->AddVector(min);
    817     Translate(min);
    818     Center.Zero();
    819   }
    820   delete(min);
    821 //  *out << Verbose(3) << "End of CenterEdge." << endl;
    822 };
    823 
    824 /** Centers the center of the atoms at (0,0,0).
    825  * \param *out output stream for debugging
    826  * \param *center return vector for translation vector
    827  */
    828 void molecule::CenterOrigin(ofstream *out)
    829 {
    830   int Num = 0;
    831   atom *ptr = start->next;  // start at first in list
    832 
    833   Center.Zero();
    834 
    835   if (ptr != end) {   //list not empty?
    836     while (ptr->next != end) {  // continue with second if present
    837       ptr = ptr->next;
    838       Num++;
    839       Center.AddVector(&ptr->x);
    840     }
    841     Center.Scale(-1./Num); // divide through total number (and sign for direction)
    842     Translate(&Center);
    843     Center.Zero();
    844   }
    845 };
    846 
    847 /** Returns vector pointing to center of all atoms.
    848  * \param *out output stream for debugging
    849  * \return pointer to center of all vector
    850  */
    851 Vector * molecule::DetermineCenterOfAll(ofstream *out)
    852 {
    853   atom *ptr = start->next;  // start at first in list
    854   Vector *a = new Vector();
    855   Vector tmp;
    856   double Num = 0;
    857 
    858   a->Zero();
    859 
    860   if (ptr != end) {   //list not empty?
    861     while (ptr->next != end) {  // continue with second if present
    862       ptr = ptr->next;
    863       Num += 1.;
    864       tmp.CopyVector(&ptr->x);
    865       a->AddVector(&tmp);
    866     }
    867     a->Scale(1./Num); // divide through total mass (and sign for direction)
    868   }
    869   //cout << Verbose(1) << "Resulting center of gravity: ";
    870   //a->Output(out);
    871   //cout << endl;
    872   return a;
    873 };
    874 
    875 /** Returns vector pointing to center of gravity.
    876  * \param *out output stream for debugging
    877  * \return pointer to center of gravity vector
    878  */
    879 Vector * molecule::DetermineCenterOfGravity(ofstream *out)
    880 {
    881   atom *ptr = start->next;  // start at first in list
    882   Vector *a = new Vector();
    883   Vector tmp;
    884   double Num = 0;
    885 
    886   a->Zero();
    887 
    888   if (ptr != end) {   //list not empty?
    889     while (ptr->next != end) {  // continue with second if present
    890       ptr = ptr->next;
    891       Num += ptr->type->mass;
    892       tmp.CopyVector(&ptr->x);
    893       tmp.Scale(ptr->type->mass);  // scale by mass
    894       a->AddVector(&tmp);
    895     }
    896     a->Scale(-1./Num); // divide through total mass (and sign for direction)
    897   }
    898 //  *out << Verbose(1) << "Resulting center of gravity: ";
    899 //  a->Output(out);
    900 //  *out << endl;
    901   return a;
    902 };
    903 
    904 /** Centers the center of gravity of the atoms at (0,0,0).
    905  * \param *out output stream for debugging
    906  * \param *center return vector for translation vector
    907  */
    908 void molecule::CenterPeriodic(ofstream *out)
    909 {
    910   DeterminePeriodicCenter(Center);
    911 };
    912 
    913 
    914 /** Centers the center of gravity of the atoms at (0,0,0).
    915  * \param *out output stream for debugging
    916  * \param *center return vector for translation vector
    917  */
    918 void molecule::CenterAtVector(ofstream *out, Vector *newcenter)
    919 {
    920   Center.CopyVector(newcenter);
    921 };
    922 
    923 
    924 /** Scales all atoms by \a *factor.
    925  * \param *factor pointer to scaling factor
    926  */
    927 void molecule::Scale(double **factor)
    928 {
    929   atom *ptr = start;
    930 
    931   while (ptr->next != end) {
    932     ptr = ptr->next;
    933     for (int j=0;j<MDSteps;j++)
    934       Trajectories[ptr].R.at(j).Scale(factor);
    935     ptr->x.Scale(factor);
    936   }
    937 };
    938 
    939 /** Translate all atoms by given vector.
    940  * \param trans[] translation vector.
    941  */
    942 void molecule::Translate(const Vector *trans)
    943 {
    944   atom *ptr = start;
    945 
    946   while (ptr->next != end) {
    947     ptr = ptr->next;
    948     for (int j=0;j<MDSteps;j++)
    949       Trajectories[ptr].R.at(j).Translate(trans);
    950     ptr->x.Translate(trans);
    951   }
    952 };
    953 
    954 /** Translate the molecule periodically in the box.
    955  * \param trans[] translation vector.
    956  */
    957 void molecule::TranslatePeriodically(const Vector *trans)
    958 {
    959   atom *ptr = NULL;
    960   Vector x;
    961   double *M = ReturnFullMatrixforSymmetric(cell_size);
    962   double *Minv = x.InverseMatrix(M);
    963   double value;
    964 
    965   // go through all atoms
    966   ptr = start->next;  // start at first in list
    967   while (ptr != end) {  // continue with second if present
    968     //ptr->Output(1,1,out);
    969     // multiply its vector with matrix inverse
    970     x.CopyVector(&ptr->x);
    971     x.Translate(trans);
    972     x.MatrixMultiplication(Minv);
    973     // truncate to [0,1] for each axis
    974     for (int i=0;i<NDIM;i++) {
    975       value = floor(fabs(x.x[i]));  // next lower integer
    976       if (x.x[i] >=0) {
    977         x.x[i] -= value;
    978       } else {
    979         x.x[i] += value+1;
    980       }
    981     }
    982     // matrix multiply
    983     x.MatrixMultiplication(M);
    984     ptr->x.CopyVector(&x);
    985     for (int j=0;j<MDSteps;j++) {
    986       x.CopyVector(&Trajectories[ptr].R.at(j));
    987       x.Translate(trans);
    988       x.MatrixMultiplication(Minv);
    989       // truncate to [0,1] for each axis
    990       for (int i=0;i<NDIM;i++) {
    991         value = floor(x.x[i]);  // next lower integer
    992         if (x.x[i] >=0) {
    993           x.x[i] -= value;
    994         } else {
    995           x.x[i] += value+1;
    996         }
    997       }
    998       // matrix multiply
    999       x.MatrixMultiplication(M);
    1000       Trajectories[ptr].R.at(j).CopyVector(&x);
    1001     }
    1002     ptr = ptr->next;
    1003   }
    1004   delete(M);
    1005   delete(Minv);
    1006 };
    1007 
    1008 
    1009 /** Mirrors all atoms against a given plane.
    1010  * \param n[] normal vector of mirror plane.
    1011  */
    1012 void molecule::Mirror(const Vector *n)
    1013 {
    1014   atom *ptr = start;
    1015 
    1016   while (ptr->next != end) {
    1017     ptr = ptr->next;
    1018     for (int j=0;j<MDSteps;j++)
    1019       Trajectories[ptr].R.at(j).Mirror(n);
    1020     ptr->x.Mirror(n);
    1021   }
    1022 };
    1023 
    1024 /** Determines center of molecule (yet not considering atom masses).
    1025  * \param center reference to return vector
    1026  */
    1027 void molecule::DeterminePeriodicCenter(Vector &center)
    1028 {
    1029   atom *Walker = start;
    1030   bond *Binder = NULL;
    1031   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    1032   double tmp;
    1033   bool flag;
    1034   Vector Testvector, Translationvector;
    1035 
    1036   do {
    1037     Center.Zero();
    1038     flag = true;
    1039     while (Walker->next != end) {
    1040       Walker = Walker->next;
    1041 #ifdef ADDHYDROGEN
    1042       if (Walker->type->Z != 1) {
    1043 #endif
    1044         Testvector.CopyVector(&Walker->x);
    1045         Testvector.InverseMatrixMultiplication(matrix);
    1046         Translationvector.Zero();
    1047         for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
    1048           Binder = ListOfBondsPerAtom[Walker->nr][i];
    1049           if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    1050             for (int j=0;j<NDIM;j++) {
    1051               tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
    1052               if ((fabs(tmp)) > BondDistance) {
    1053                 flag = false;
    1054                 cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl;
    1055                 if (tmp > 0)
    1056                   Translationvector.x[j] -= 1.;
    1057                 else
    1058                   Translationvector.x[j] += 1.;
    1059               }
    1060             }
    1061         }
    1062         Testvector.AddVector(&Translationvector);
    1063         Testvector.MatrixMultiplication(matrix);
    1064         Center.AddVector(&Testvector);
    1065         cout << Verbose(1) << "vector is: ";
    1066         Testvector.Output((ofstream *)&cout);
    1067         cout << endl;
    1068 #ifdef ADDHYDROGEN
    1069         // now also change all hydrogens
    1070         for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
    1071           Binder = ListOfBondsPerAtom[Walker->nr][i];
    1072           if (Binder->GetOtherAtom(Walker)->type->Z == 1) {
    1073             Testvector.CopyVector(&Binder->GetOtherAtom(Walker)->x);
    1074             Testvector.InverseMatrixMultiplication(matrix);
    1075             Testvector.AddVector(&Translationvector);
    1076             Testvector.MatrixMultiplication(matrix);
    1077             Center.AddVector(&Testvector);
    1078             cout << Verbose(1) << "Hydrogen vector is: ";
    1079             Testvector.Output((ofstream *)&cout);
    1080             cout << endl;
    1081           }
    1082         }
    1083       }
    1084 #endif
    1085     }
    1086   } while (!flag);
    1087   Free(&matrix);
    1088   Center.Scale(1./(double)AtomCount);
    1089 };
    1090 
    1091 /** Transforms/Rotates the given molecule into its principal axis system.
    1092  * \param *out output stream for debugging
    1093  * \param DoRotate whether to rotate (true) or only to determine the PAS.
    1094  */
    1095 void molecule::PrincipalAxisSystem(ofstream *out, bool DoRotate)
    1096 {
    1097   atom *ptr = start;  // start at first in list
    1098   double InertiaTensor[NDIM*NDIM];
    1099   Vector *CenterOfGravity = DetermineCenterOfGravity(out);
    1100 
    1101   CenterPeriodic(out);
    1102 
    1103   // reset inertia tensor
    1104   for(int i=0;i<NDIM*NDIM;i++)
    1105     InertiaTensor[i] = 0.;
    1106 
    1107   // sum up inertia tensor
    1108   while (ptr->next != end) {
    1109     ptr = ptr->next;
    1110     Vector x;
    1111     x.CopyVector(&ptr->x);
    1112     //x.SubtractVector(CenterOfGravity);
    1113     InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
    1114     InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
    1115     InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
    1116     InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
    1117     InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
    1118     InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
    1119     InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
    1120     InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
    1121     InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
    1122   }
    1123   // print InertiaTensor for debugging
    1124   *out << "The inertia tensor is:" << endl;
    1125   for(int i=0;i<NDIM;i++) {
    1126     for(int j=0;j<NDIM;j++)
    1127       *out << InertiaTensor[i*NDIM+j] << " ";
    1128     *out << endl;
    1129   }
    1130   *out << endl;
    1131 
    1132   // diagonalize to determine principal axis system
    1133   gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
    1134   gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
    1135   gsl_vector *eval = gsl_vector_alloc(NDIM);
    1136   gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
    1137   gsl_eigen_symmv(&m.matrix, eval, evec, T);
    1138   gsl_eigen_symmv_free(T);
    1139   gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
    1140 
    1141   for(int i=0;i<NDIM;i++) {
    1142     *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
    1143     *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
    1144   }
    1145 
    1146   // check whether we rotate or not
    1147   if (DoRotate) {
    1148     *out << Verbose(1) << "Transforming molecule into PAS ... ";
    1149     // the eigenvectors specify the transformation matrix
    1150     ptr = start;
    1151     while (ptr->next != end) {
    1152       ptr = ptr->next;
    1153       for (int j=0;j<MDSteps;j++)
    1154         Trajectories[ptr].R.at(j).MatrixMultiplication(evec->data);
    1155       ptr->x.MatrixMultiplication(evec->data);
    1156     }
    1157     *out << "done." << endl;
    1158 
    1159     // summing anew for debugging (resulting matrix has to be diagonal!)
    1160     // reset inertia tensor
    1161     for(int i=0;i<NDIM*NDIM;i++)
    1162       InertiaTensor[i] = 0.;
    1163 
    1164     // sum up inertia tensor
    1165     ptr = start;
    1166     while (ptr->next != end) {
    1167       ptr = ptr->next;
    1168       Vector x;
    1169       x.CopyVector(&ptr->x);
    1170       //x.SubtractVector(CenterOfGravity);
    1171       InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
    1172       InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
    1173       InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
    1174       InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
    1175       InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
    1176       InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
    1177       InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
    1178       InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
    1179       InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
    1180     }
    1181     // print InertiaTensor for debugging
    1182     *out << "The inertia tensor is:" << endl;
    1183     for(int i=0;i<NDIM;i++) {
    1184       for(int j=0;j<NDIM;j++)
    1185         *out << InertiaTensor[i*NDIM+j] << " ";
    1186       *out << endl;
    1187     }
    1188     *out << endl;
    1189   }
    1190 
    1191   // free everything
    1192   delete(CenterOfGravity);
    1193   gsl_vector_free(eval);
    1194   gsl_matrix_free(evec);
    1195 };
    1196 
    1197 /** Evaluates the potential energy used for constrained molecular dynamics.
    1198  * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
    1199  *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}\f$ is minimum distance between
    1200  *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
    1201  * Note that for the second term we have to solve the following linear system:
    1202  * \f$-c_1 \cdot n_1 + c_2 \cdot n_2 + C \cdot n_3 = - p_2 + p_1\f$, where \f$c_1\f$, \f$c_2\f$ and \f$C\f$ are constants,
    1203  * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
    1204  * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
    1205  * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
    1206  * \param *out output stream for debugging
    1207  * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
    1208  * \param startstep start configuration (MDStep in molecule::trajectories)
    1209  * \param endstep end configuration (MDStep in molecule::trajectories)
    1210  * \param *constants constant in front of each term
    1211  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
    1212  * \return potential energy
    1213  * \note This routine is scaling quadratically which is not optimal.
    1214  * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
    1215  */
    1216 double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
    1217 {
    1218   double result = 0., tmp, Norm1, Norm2;
    1219   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
    1220   Vector trajectory1, trajectory2, normal, TestVector;
    1221   gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    1222   gsl_vector *x = gsl_vector_alloc(NDIM);
    1223 
    1224   // go through every atom
    1225   Walker = start;
    1226   while (Walker->next != end) {
    1227     Walker = Walker->next;
    1228     // first term: distance to target
    1229     Runner = PermutationMap[Walker->nr];   // find target point
    1230     tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
    1231     tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    1232     result += constants[0] * tmp;
    1233     //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
    1234    
    1235     // second term: sum of distances to other trajectories
    1236     Runner = start;
    1237     while (Runner->next != end) {
    1238       Runner = Runner->next;
    1239       if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
    1240         break;
    1241       // determine normalized trajectories direction vector (n1, n2)
    1242       Sprinter = PermutationMap[Walker->nr];   // find first target point
    1243       trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
    1244       trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
    1245       trajectory1.Normalize();
    1246       Norm1 = trajectory1.Norm();
    1247       Sprinter = PermutationMap[Runner->nr];   // find second target point
    1248       trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
    1249       trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
    1250       trajectory2.Normalize();
    1251       Norm2 = trajectory1.Norm();
    1252       // check whether either is zero()
    1253       if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    1254         tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
    1255       } else if (Norm1 < MYEPSILON) {
    1256         Sprinter = PermutationMap[Walker->nr];   // find first target point
    1257         trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
    1258         trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
    1259         trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    1260         trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
    1261         tmp = trajectory1.Norm();  // remaining norm is distance
    1262       } else if (Norm2 < MYEPSILON) {
    1263         Sprinter = PermutationMap[Runner->nr];   // find second target point
    1264         trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
    1265         trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
    1266         trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    1267         trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
    1268         tmp = trajectory2.Norm();  // remaining norm is distance
    1269       } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
    1270 //        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
    1271 //        *out << trajectory1;
    1272 //        *out << " and ";
    1273 //        *out << trajectory2;
    1274         tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
    1275 //        *out << " with distance " << tmp << "." << endl;
    1276       } else { // determine distance by finding minimum distance
    1277 //        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
    1278 //        *out << endl;
    1279 //        *out << "First Trajectory: ";
    1280 //        *out << trajectory1 << endl;
    1281 //        *out << "Second Trajectory: ";
    1282 //        *out << trajectory2 << endl;
    1283         // determine normal vector for both
    1284         normal.MakeNormalVector(&trajectory1, &trajectory2);
    1285         // print all vectors for debugging
    1286 //        *out << "Normal vector in between: ";
    1287 //        *out << normal << endl;
    1288         // setup matrix
    1289         for (int i=NDIM;i--;) {
    1290           gsl_matrix_set(A, 0, i, trajectory1.x[i]);
    1291           gsl_matrix_set(A, 1, i, trajectory2.x[i]);
    1292           gsl_matrix_set(A, 2, i, normal.x[i]);
    1293           gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
    1294         }
    1295         // solve the linear system by Householder transformations
    1296         gsl_linalg_HH_svx(A, x);
    1297         // distance from last component
    1298         tmp = gsl_vector_get(x,2);
    1299 //        *out << " with distance " << tmp << "." << endl;
    1300         // test whether we really have the intersection (by checking on c_1 and c_2)
    1301         TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
    1302         trajectory2.Scale(gsl_vector_get(x,1));
    1303         TestVector.AddVector(&trajectory2);
    1304         normal.Scale(gsl_vector_get(x,2));
    1305         TestVector.AddVector(&normal);
    1306         TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
    1307         trajectory1.Scale(gsl_vector_get(x,0));
    1308         TestVector.SubtractVector(&trajectory1);
    1309         if (TestVector.Norm() < MYEPSILON) {
    1310 //          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl; 
    1311         } else {
    1312 //          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
    1313 //          *out << TestVector;
    1314 //          *out << "." << endl; 
    1315         }
    1316       }
    1317       // add up
    1318       tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    1319       if (fabs(tmp) > MYEPSILON) {
    1320         result += constants[1] * 1./tmp;
    1321         //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
    1322       }
    1323     }
    1324      
    1325     // third term: penalty for equal targets
    1326     Runner = start;
    1327     while (Runner->next != end) {
    1328       Runner = Runner->next;
    1329       if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
    1330         Sprinter = PermutationMap[Walker->nr];
    1331 //        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
    1332 //        *out << Trajectories[Sprinter].R.at(endstep);
    1333 //        *out << ", penalting." << endl;
    1334         result += constants[2];
    1335         //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
    1336       }
    1337     }
    1338   }
    1339  
    1340   return result;
    1341 };
    1342 
    1343 void PrintPermutationMap(ofstream *out, atom **PermutationMap, int Nr)
    1344 {
    1345   stringstream zeile1, zeile2;
    1346   int *DoubleList = Malloc<int>(Nr, "PrintPermutationMap: *DoubleList");
    1347   int doubles = 0;
    1348   for (int i=0;i<Nr;i++)
    1349     DoubleList[i] = 0;
    1350   zeile1 << "PermutationMap: ";
    1351   zeile2 << "                ";
    1352   for (int i=0;i<Nr;i++) {
    1353     DoubleList[PermutationMap[i]->nr]++;
    1354     zeile1 << i << " ";
    1355     zeile2 << PermutationMap[i]->nr << " ";
    1356   }
    1357   for (int i=0;i<Nr;i++)
    1358     if (DoubleList[i] > 1)
    1359     doubles++;
    1360  // *out << "Found " << doubles << " Doubles." << endl;
    1361   Free(&DoubleList);
    1362 //  *out << zeile1.str() << endl << zeile2.str() << endl;
    1363 };
    1364 
    1365 /** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
    1366  * We do the following:
    1367  *  -# Generate a distance list from all source to all target points
    1368  *  -# Sort this per source point
    1369  *  -# Take for each source point the target point with minimum distance, use this as initial permutation
    1370  *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
    1371  *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
    1372  *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
    1373  *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
    1374  *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
    1375  *     if this decreases the conditional potential.
    1376  *  -# finished.
    1377  *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
    1378  *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
    1379  *     right direction).
    1380  *  -# Finally, we calculate the potential energy and return.
    1381  * \param *out output stream for debugging
    1382  * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
    1383  * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
    1384  * \param endstep step giving final position in constrained MD
    1385  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
    1386  * \sa molecule::VerletForceIntegration()
    1387  * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
    1388  * \todo The constrained potential's constants are set to fixed values right now, but they should scale based on checks of the system in order
    1389  *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
    1390  * \bug this all is not O(N log N) but O(N^2)
    1391  */
    1392 double molecule::MinimiseConstrainedPotential(ofstream *out, atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
    1393 {
    1394   double Potential, OldPotential, OlderPotential;
    1395   PermutationMap = Malloc<atom*>(AtomCount, "molecule::MinimiseConstrainedPotential: **PermutationMap");
    1396   DistanceMap **DistanceList = Malloc<DistanceMap*>(AtomCount, "molecule::MinimiseConstrainedPotential: **DistanceList");
    1397   DistanceMap::iterator *DistanceIterators = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
    1398   int *DoubleList = Malloc<int>(AtomCount, "molecule::MinimiseConstrainedPotential: *DoubleList");
    1399   DistanceMap::iterator *StepList = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: *StepList");
    1400   double constants[3];
    1401   int round;
    1402   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
    1403   DistanceMap::iterator Rider, Strider;
    1404  
    1405   /// Minimise the potential
    1406   // set Lagrange multiplier constants
    1407   constants[0] = 10.;
    1408   constants[1] = 1.;
    1409   constants[2] = 1e+7;    // just a huge penalty
    1410   // generate the distance list
    1411   *out << Verbose(1) << "Creating the distance list ... " << endl;
    1412   for (int i=AtomCount; i--;) {
    1413     DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
    1414     DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    1415     DistanceList[i]->clear();
    1416   }
    1417   *out << Verbose(1) << "Filling the distance list ... " << endl;
    1418   Walker = start;
    1419   while (Walker->next != end) {
    1420     Walker = Walker->next;
    1421     Runner = start;
    1422     while(Runner->next != end) {
    1423       Runner = Runner->next;
    1424       DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
    1425     }
    1426   }
    1427   // create the initial PermutationMap (source -> target)
    1428   Walker = start;
    1429   while (Walker->next != end) {
    1430     Walker = Walker->next;
    1431     StepList[Walker->nr] = DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
    1432     PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
    1433     DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
    1434     DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
    1435     *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
    1436   }
    1437   *out << Verbose(1) << "done." << endl;
    1438   // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
    1439   *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
    1440   Walker = start;
    1441   DistanceMap::iterator NewBase;
    1442   OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
    1443   while ((OldPotential) > constants[2]) {
    1444     PrintPermutationMap(out, PermutationMap, AtomCount);
    1445     Walker = Walker->next;
    1446     if (Walker == end) // round-robin at the end
    1447       Walker = start->next;
    1448     if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
    1449       continue;
    1450     // now, try finding a new one
    1451     NewBase = DistanceIterators[Walker->nr];  // store old base
    1452     do {
    1453       NewBase++;  // take next further distance in distance to targets list that's a target of no one
    1454     } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
    1455     if (NewBase != DistanceList[Walker->nr]->end()) {
    1456       PermutationMap[Walker->nr] = NewBase->second;
    1457       Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
    1458       if (Potential > OldPotential) { // undo
    1459         PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
    1460       } else {  // do
    1461         DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
    1462         DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
    1463         DistanceIterators[Walker->nr] = NewBase;
    1464         OldPotential = Potential;
    1465         *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
    1466       }
    1467     }
    1468   }
    1469   for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
    1470     if (DoubleList[i] > 1) {
    1471       cerr << "Failed to create an injective PermutationMap!" << endl;
    1472       exit(1);
    1473     }
    1474   *out << Verbose(1) << "done." << endl;
    1475   Free(&DoubleList);
    1476   // argument minimise the constrained potential in this injective PermutationMap
    1477   *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
    1478   OldPotential = 1e+10;
    1479   round = 0;
    1480   do {
    1481     *out << "Starting round " << ++round << " ... " << endl;
    1482     OlderPotential = OldPotential;
    1483     do {
    1484       Walker = start;
    1485       while (Walker->next != end) { // pick one
    1486         Walker = Walker->next;
    1487         PrintPermutationMap(out, PermutationMap, AtomCount);
    1488         Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
    1489         Strider = DistanceIterators[Walker->nr];  //remember old iterator
    1490         DistanceIterators[Walker->nr] = StepList[Walker->nr]; 
    1491         if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
    1492           DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->begin();
    1493           break;
    1494         }
    1495         //*out << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
    1496         // find source of the new target
    1497         Runner = start->next;
    1498         while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    1499           if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
    1500             //*out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
    1501             break;
    1502           }
    1503           Runner = Runner->next;
    1504         }
    1505         if (Runner != end) { // we found the other source
    1506           // then look in its distance list for Sprinter
    1507           Rider = DistanceList[Runner->nr]->begin();
    1508           for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
    1509             if (Rider->second == Sprinter)
    1510               break;
    1511           if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
    1512             //*out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
    1513             // exchange both
    1514             PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
    1515             PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
    1516             PrintPermutationMap(out, PermutationMap, AtomCount);
    1517             // calculate the new potential
    1518             //*out << Verbose(2) << "Checking new potential ..." << endl;
    1519             Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
    1520             if (Potential > OldPotential) { // we made everything worse! Undo ...
    1521               //*out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    1522               //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *DistanceIterators[Runner->nr]->second << "." << endl;
    1523               // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    1524               PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
    1525               // Undo for Walker
    1526               DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
    1527               //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *DistanceIterators[Walker->nr]->second << "." << endl;
    1528               PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
    1529             } else {
    1530               DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
    1531               *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
    1532               OldPotential = Potential;
    1533             }
    1534             if (Potential > constants[2]) {
    1535               cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
    1536               exit(255);
    1537             }
    1538             //*out << endl;
    1539           } else {
    1540             cerr << "ERROR: " << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
    1541             exit(255);
    1542           }
    1543         } else {
    1544           PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // new target has no source!
    1545         }
    1546         StepList[Walker->nr]++; // take next farther distance target
    1547       }
    1548     } while (Walker->next != end);
    1549   } while ((OlderPotential - OldPotential) > 1e-3);
    1550   *out << Verbose(1) << "done." << endl;
    1551 
    1552  
    1553   /// free memory and return with evaluated potential
    1554   for (int i=AtomCount; i--;)
    1555     DistanceList[i]->clear();
    1556   Free(&DistanceList);
    1557   Free(&DistanceIterators);
    1558   return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
    1559 };
    1560 
    1561 /** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
    1562  * \param *out output stream for debugging
    1563  * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
    1564  * \param endstep step giving final position in constrained MD
    1565  * \param **PermutationMap mapping between the atom label of the initial and the final configuration
    1566  * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
    1567  * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
    1568  */
    1569 void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
    1570 {
    1571   double constant = 10.;
    1572   atom *Walker = NULL, *Sprinter = NULL;
    1573  
    1574   /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
    1575   *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
    1576   Walker = start;
    1577   while (Walker->next != NULL) {
    1578     Walker = Walker->next;
    1579     Sprinter = PermutationMap[Walker->nr];
    1580     // set forces
    1581     for (int i=NDIM;i++;)
    1582       Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
    1583   } 
    1584   *out << Verbose(1) << "done." << endl;
    1585 };
    1586 
    1587 /** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
    1588  * Note, step number is config::MaxOuterStep
    1589  * \param *out output stream for debugging
    1590  * \param startstep stating initial configuration in molecule::Trajectories
    1591  * \param endstep stating final configuration in molecule::Trajectories
    1592  * \param &config configuration structure
    1593  * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
    1594  * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
    1595  */
    1596 bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
    1597 {
    1598   molecule *mol = NULL;
    1599   bool status = true;
    1600   int MaxSteps = configuration.MaxOuterStep;
    1601   MoleculeListClass *MoleculePerStep = new MoleculeListClass();
    1602   // Get the Permutation Map by MinimiseConstrainedPotential
    1603   atom **PermutationMap = NULL;
    1604   atom *Walker = NULL, *Sprinter = NULL;
    1605   if (!MapByIdentity)
    1606     MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
    1607   else {
    1608     PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
    1609     Walker = start;
    1610     while (Walker->next != end) {
    1611       Walker = Walker->next;
    1612       PermutationMap[Walker->nr] = Walker;   // always pick target with the smallest distance
    1613     }
    1614   }
    1615 
    1616   // check whether we have sufficient space in Trajectories for each atom
    1617   Walker = start;
    1618   while (Walker->next != end) {
    1619     Walker = Walker->next;
    1620     if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
    1621       //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
    1622       Trajectories[Walker].R.resize(MaxSteps+1);
    1623       Trajectories[Walker].U.resize(MaxSteps+1);
    1624       Trajectories[Walker].F.resize(MaxSteps+1);
    1625     }
    1626   }
    1627   // push endstep to last one
    1628   Walker = start;
    1629   while (Walker->next != end) { // remove the endstep (is now the very last one)
    1630     Walker = Walker->next;
    1631     for (int n=NDIM;n--;) {
    1632       Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
    1633       Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
    1634       Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
    1635     }
    1636   }
    1637   endstep = MaxSteps;
    1638  
    1639   // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
    1640   *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
    1641   for (int step = 0; step <= MaxSteps; step++) {
    1642     mol = new molecule(elemente);
    1643     MoleculePerStep->insert(mol);
    1644     Walker = start;
    1645     while (Walker->next != end) {
    1646       Walker = Walker->next;
    1647       // add to molecule list
    1648       Sprinter = mol->AddCopyAtom(Walker);
    1649       for (int n=NDIM;n--;) {
    1650         Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
    1651         // add to Trajectories
    1652         //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    1653         if (step < MaxSteps) {
    1654           Trajectories[Walker].R.at(step).x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
    1655           Trajectories[Walker].U.at(step).x[n] = 0.;
    1656           Trajectories[Walker].F.at(step).x[n] = 0.;
    1657         }
    1658       }
    1659     }
    1660   }
    1661   MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit
    1662  
    1663   // store the list to single step files
    1664   int *SortIndex = Malloc<int>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
    1665   for (int i=AtomCount; i--; )
    1666     SortIndex[i] = i;
    1667   status = MoleculePerStep->OutputConfigForListOfFragments(out, &configuration, SortIndex);
    1668  
    1669   // free and return
    1670   Free(&PermutationMap);
    1671   delete(MoleculePerStep);
    1672   return status;
    1673 };
    1674 
    1675 /** Parses nuclear forces from file and performs Verlet integration.
    1676  * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
    1677  * have to transform them).
    1678  * This adds a new MD step to the config file.
    1679  * \param *out output stream for debugging
    1680  * \param *file filename
    1681  * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
    1682  * \param delta_t time step width in atomic units
    1683  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
    1684  * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential()
    1685  * \return true - file found and parsed, false - file not found or imparsable
    1686  * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
    1687  */
    1688 bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
    1689 {
    1690   atom *walker = NULL;
    1691   ifstream input(file);
    1692   string token;
    1693   stringstream item;
    1694   double IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
    1695   ForceMatrix Force;
    1696 
    1697   CountElements();  // make sure ElementsInMolecule is up to date
    1698 
    1699   // check file
    1700   if (input == NULL) {
    1701     return false;
    1702   } else {
    1703     // parse file into ForceMatrix
    1704     if (!Force.ParseMatrix(file, 0,0,0)) {
    1705       cerr << "Could not parse Force Matrix file " << file << "." << endl;
    1706       return false;
    1707     }
    1708     if (Force.RowCounter[0] != AtomCount) {
    1709       cerr << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl;
    1710       return false;
    1711     }
    1712     // correct Forces
    1713     for(int d=0;d<NDIM;d++)
    1714       Vector[d] = 0.;
    1715     for(int i=0;i<AtomCount;i++)
    1716       for(int d=0;d<NDIM;d++) {
    1717         Vector[d] += Force.Matrix[0][i][d+5];
    1718       }
    1719     for(int i=0;i<AtomCount;i++)
    1720       for(int d=0;d<NDIM;d++) {
    1721         Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
    1722       }
    1723     // solve a constrained potential if we are meant to
    1724     if (configuration.DoConstrainedMD) {
    1725       // calculate forces and potential
    1726       atom **PermutationMap = NULL;
    1727       ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    1728       EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
    1729       Free(&PermutationMap);
    1730     }
    1731    
    1732     // and perform Verlet integration for each atom with position, velocity and force vector
    1733     walker = start;
    1734     while (walker->next != end) { // go through every atom of this element
    1735       walker = walker->next;
    1736       //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    1737       // check size of vectors
    1738       if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
    1739         //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
    1740         Trajectories[walker].R.resize(MDSteps+10);
    1741         Trajectories[walker].U.resize(MDSteps+10);
    1742         Trajectories[walker].F.resize(MDSteps+10);
    1743       }
    1744      
    1745       // Update R (and F)
    1746       for (int d=0; d<NDIM; d++) {
    1747         Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
    1748         Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
    1749         Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
    1750         Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
    1751       }
    1752       // Update U
    1753       for (int d=0; d<NDIM; d++) {
    1754         Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
    1755         Trajectories[walker].U.at(MDSteps).x[d] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
    1756       }
    1757 //      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
    1758 //      for (int d=0;d<NDIM;d++)
    1759 //        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
    1760 //      out << ")\t(";
    1761 //      for (int d=0;d<NDIM;d++)
    1762 //        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
    1763 //      out << ")" << endl;
    1764             // next atom
    1765     }
    1766   }
    1767   // correct velocities (rather momenta) so that center of mass remains motionless
    1768   for(int d=0;d<NDIM;d++)
    1769     Vector[d] = 0.;
    1770   IonMass = 0.;
    1771   walker = start;
    1772   while (walker->next != end) { // go through every atom
    1773     walker = walker->next;
    1774     IonMass += walker->type->mass;  // sum up total mass
    1775     for(int d=0;d<NDIM;d++) {
    1776       Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
    1777     }
    1778   }
    1779   // correct velocities (rather momenta) so that center of mass remains motionless
    1780   for(int d=0;d<NDIM;d++)
    1781     Vector[d] /= IonMass;
    1782   ActualTemp = 0.;
    1783   walker = start;
    1784   while (walker->next != end) { // go through every atom of this element
    1785     walker = walker->next;
    1786     for(int d=0;d<NDIM;d++) {
    1787       Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
    1788       ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
    1789     }
    1790   }
    1791   Thermostats(configuration, ActualTemp, Berendsen);
    1792   MDSteps++;
    1793 
    1794 
    1795   // exit
    1796   return true;
    1797 };
    1798 
    1799 /** Implementation of various thermostats.
    1800  * All these thermostats apply an additional force which has the following forms:
    1801  * -# Woodcock
    1802  *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
    1803  * -# Gaussian
    1804  *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
    1805  * -# Langevin
    1806  *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$ 
    1807  * -# Berendsen
    1808  *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
    1809  * -# Nose-Hoover
    1810  *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
    1811  * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
    1812  * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
    1813  * belatedly and the constraint force set.
    1814  * \param *P Problem at hand
    1815  * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
    1816  * \sa InitThermostat()
    1817  */
    1818 void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
    1819 {
    1820   double ekin = 0.;
    1821   double E = 0., G = 0.;
    1822   double delta_alpha = 0.;
    1823   double ScaleTempFactor;
    1824   double sigma;
    1825   double IonMass;
    1826   int d;
    1827   gsl_rng * r;
    1828   const gsl_rng_type * T;
    1829   double *U = NULL, *F = NULL, FConstraint[NDIM];
    1830   atom *walker = NULL;
    1831  
    1832   // calculate scale configuration
    1833   ScaleTempFactor = configuration.TargetTemp/ActualTemp;
    1834    
    1835   // differentating between the various thermostats
    1836   switch(Thermostat) {
    1837      case None:
    1838       cout << Verbose(2) <<  "Applying no thermostat..." << endl;
    1839       break;
    1840      case Woodcock:
    1841       if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
    1842         cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
    1843         walker = start;
    1844         while (walker->next != end) { // go through every atom of this element
    1845           walker = walker->next;
    1846           IonMass = walker->type->mass;
    1847           U = Trajectories[walker].U.at(MDSteps).x;
    1848           if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
    1849             for (d=0; d<NDIM; d++) {
    1850               U[d] *= sqrt(ScaleTempFactor);
    1851               ekin += 0.5*IonMass * U[d]*U[d];
    1852             }
    1853         }
    1854       }
    1855       break;
    1856      case Gaussian:
    1857       cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
    1858       walker = start;
    1859       while (walker->next != end) { // go through every atom of this element
    1860         walker = walker->next;
    1861         IonMass = walker->type->mass;
    1862         U = Trajectories[walker].U.at(MDSteps).x;
    1863         F = Trajectories[walker].F.at(MDSteps).x;
    1864         if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
    1865           for (d=0; d<NDIM; d++) {
    1866             G += U[d] * F[d];
    1867             E += U[d]*U[d]*IonMass;
    1868           }
    1869       }
    1870       cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
    1871       walker = start;
    1872       while (walker->next != end) { // go through every atom of this element
    1873         walker = walker->next;
    1874         IonMass = walker->type->mass;
    1875         U = Trajectories[walker].U.at(MDSteps).x;
    1876         F = Trajectories[walker].F.at(MDSteps).x;
    1877         if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
    1878           for (d=0; d<NDIM; d++) {
    1879             FConstraint[d] = (G/E) * (U[d]*IonMass);
    1880             U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
    1881             ekin += IonMass * U[d]*U[d];
    1882           }
    1883       }
    1884       break;
    1885      case Langevin:
    1886       cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
    1887       // init random number generator
    1888       gsl_rng_env_setup();
    1889       T = gsl_rng_default;
    1890       r = gsl_rng_alloc (T);
    1891       // Go through each ion
    1892       walker = start;
    1893       while (walker->next != end) { // go through every atom of this element
    1894         walker = walker->next;
    1895         IonMass = walker->type->mass;
    1896         sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
    1897         U = Trajectories[walker].U.at(MDSteps).x;
    1898         F = Trajectories[walker].F.at(MDSteps).x;
    1899         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    1900           // throw a dice to determine whether it gets hit by a heat bath particle
    1901           if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) { 
    1902             cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
    1903             // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
    1904             for (d=0; d<NDIM; d++) {
    1905               U[d] = gsl_ran_gaussian (r, sigma);
    1906             }
    1907             cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
    1908           }
    1909           for (d=0; d<NDIM; d++)
    1910             ekin += 0.5*IonMass * U[d]*U[d];
    1911         }
    1912       }
    1913       break;
    1914      case Berendsen:
    1915       cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
    1916       walker = start;
    1917       while (walker->next != end) { // go through every atom of this element
    1918         walker = walker->next;
    1919         IonMass = walker->type->mass;
    1920         U = Trajectories[walker].U.at(MDSteps).x;
    1921         F = Trajectories[walker].F.at(MDSteps).x;
    1922         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    1923           for (d=0; d<NDIM; d++) {
    1924             U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
    1925             ekin += 0.5*IonMass * U[d]*U[d];
    1926           }
    1927         }
    1928       }
    1929       break;
    1930      case NoseHoover:
    1931       cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
    1932       // dynamically evolve alpha (the additional degree of freedom)
    1933       delta_alpha = 0.;
    1934       walker = start;
    1935       while (walker->next != end) { // go through every atom of this element
    1936         walker = walker->next;
    1937         IonMass = walker->type->mass;
    1938         U = Trajectories[walker].U.at(MDSteps).x;
    1939         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    1940           for (d=0; d<NDIM; d++) {
    1941             delta_alpha += U[d]*U[d]*IonMass;
    1942           }
    1943         }
    1944       }
    1945       delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
    1946       configuration.alpha += delta_alpha*configuration.Deltat;
    1947       cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
    1948       // apply updated alpha as additional force
    1949       walker = start;
    1950       while (walker->next != end) { // go through every atom of this element
    1951         walker = walker->next;
    1952         IonMass = walker->type->mass;
    1953         U = Trajectories[walker].U.at(MDSteps).x;
    1954         if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
    1955           for (d=0; d<NDIM; d++) {
    1956               FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
    1957               U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
    1958               ekin += (0.5*IonMass) * U[d]*U[d];
    1959             }
    1960         }
    1961       }
    1962       break;
    1963   }   
    1964   cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
    1965 };
    1966 
    1967 /** Align all atoms in such a manner that given vector \a *n is along z axis.
    1968  * \param n[] alignment vector.
    1969  */
    1970 void molecule::Align(Vector *n)
    1971 {
    1972   atom *ptr = start;
    1973   double alpha, tmp;
    1974   Vector z_axis;
    1975   z_axis.x[0] = 0.;
    1976   z_axis.x[1] = 0.;
    1977   z_axis.x[2] = 1.;
    1978 
    1979   // rotate on z-x plane
    1980   cout << Verbose(0) << "Begin of Aligning all atoms." << endl;
    1981   alpha = atan(-n->x[0]/n->x[2]);
    1982   cout << Verbose(1) << "Z-X-angle: " << alpha << " ... ";
    1983   while (ptr->next != end) {
    1984     ptr = ptr->next;
    1985     tmp = ptr->x.x[0];
    1986     ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
    1987     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
    1988     for (int j=0;j<MDSteps;j++) {
    1989       tmp = Trajectories[ptr].R.at(j).x[0];
    1990       Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
    1991       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
    1992     }
    1993   }
    1994   // rotate n vector
    1995   tmp = n->x[0];
    1996   n->x[0] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
    1997   n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    1998   cout << Verbose(1) << "alignment vector after first rotation: ";
    1999   n->Output((ofstream *)&cout);
    2000   cout << endl;
    2001 
    2002   // rotate on z-y plane
    2003   ptr = start;
    2004   alpha = atan(-n->x[1]/n->x[2]);
    2005   cout << Verbose(1) << "Z-Y-angle: " << alpha << " ... ";
    2006   while (ptr->next != end) {
    2007     ptr = ptr->next;
    2008     tmp = ptr->x.x[1];
    2009     ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
    2010     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
    2011     for (int j=0;j<MDSteps;j++) {
    2012       tmp = Trajectories[ptr].R.at(j).x[1];
    2013       Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
    2014       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
    2015     }
    2016   }
    2017   // rotate n vector (for consistency check)
    2018   tmp = n->x[1];
    2019   n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
    2020   n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    2021 
    2022   cout << Verbose(1) << "alignment vector after second rotation: ";
    2023   n->Output((ofstream *)&cout);
    2024   cout << Verbose(1) << endl;
    2025   cout << Verbose(0) << "End of Aligning all atoms." << endl;
    2026 };
    2027 
    2028625/** Removes atom from molecule list and deletes it.
    2029626 * \param *pointer atom to be removed
    … …  
    2117714  //return result;
    2118715  return true; /// probably not gonna use the check no more
    2119 };
    2120 
    2121 /** Calculates sum over least square distance to line hidden in \a *x.
    2122  * \param *x offset and direction vector
    2123  * \param *params pointer to lsq_params structure
    2124  * \return \f$ sum_i^N | y_i - (a + t_i b)|^2\f$
    2125  */
    2126 double LeastSquareDistance (const gsl_vector * x, void * params)
    2127 {
    2128   double res = 0, t;
    2129   Vector a,b,c,d;
    2130   struct lsq_params *par = (struct lsq_params *)params;
    2131   atom *ptr = par->mol->start;
    2132 
    2133   // initialize vectors
    2134   a.x[0] = gsl_vector_get(x,0);
    2135   a.x[1] = gsl_vector_get(x,1);
    2136   a.x[2] = gsl_vector_get(x,2);
    2137   b.x[0] = gsl_vector_get(x,3);
    2138   b.x[1] = gsl_vector_get(x,4);
    2139   b.x[2] = gsl_vector_get(x,5);
    2140   // go through all atoms
    2141   while (ptr != par->mol->end) {
    2142     ptr = ptr->next;
    2143     if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
    2144       c.CopyVector(&ptr->x);  // copy vector to temporary one
    2145       c.SubtractVector(&a);   // subtract offset vector
    2146       t = c.ScalarProduct(&b);           // get direction parameter
    2147       d.CopyVector(&b);       // and create vector
    2148       d.Scale(&t);
    2149       c.SubtractVector(&d);   // ... yielding distance vector
    2150       res += d.ScalarProduct((const Vector *)&d);        // add squared distance
    2151     }
    2152   }
    2153   return res;
    2154 };
    2155 
    2156 /** By minimizing the least square distance gains alignment vector.
    2157  * \bug this is not yet working properly it seems
    2158  */
    2159 void molecule::GetAlignvector(struct lsq_params * par) const
    2160 {
    2161     int np = 6;
    2162 
    2163    const gsl_multimin_fminimizer_type *T =
    2164      gsl_multimin_fminimizer_nmsimplex;
    2165    gsl_multimin_fminimizer *s = NULL;
    2166    gsl_vector *ss;
    2167    gsl_multimin_function minex_func;
    2168 
    2169    size_t iter = 0, i;
    2170    int status;
    2171    double size;
    2172 
    2173    /* Initial vertex size vector */
    2174    ss = gsl_vector_alloc (np);
    2175 
    2176    /* Set all step sizes to 1 */
    2177    gsl_vector_set_all (ss, 1.0);
    2178 
    2179    /* Starting point */
    2180    par->x = gsl_vector_alloc (np);
    2181    par->mol = this;
    2182 
    2183    gsl_vector_set (par->x, 0, 0.0);  // offset
    2184    gsl_vector_set (par->x, 1, 0.0);
    2185    gsl_vector_set (par->x, 2, 0.0);
    2186    gsl_vector_set (par->x, 3, 0.0);  // direction
    2187    gsl_vector_set (par->x, 4, 0.0);
    2188    gsl_vector_set (par->x, 5, 1.0);
    2189 
    2190    /* Initialize method and iterate */
    2191    minex_func.f = &LeastSquareDistance;
    2192    minex_func.n = np;
    2193    minex_func.params = (void *)par;
    2194 
    2195    s = gsl_multimin_fminimizer_alloc (T, np);
    2196    gsl_multimin_fminimizer_set (s, &minex_func, par->x, ss);
    2197 
    2198    do
    2199      {
    2200        iter++;
    2201        status = gsl_multimin_fminimizer_iterate(s);
    2202 
    2203        if (status)
    2204          break;
    2205 
    2206        size = gsl_multimin_fminimizer_size (s);
    2207        status = gsl_multimin_test_size (size, 1e-2);
    2208 
    2209        if (status == GSL_SUCCESS)
    2210          {
    2211            printf ("converged to minimum at\n");
    2212          }
    2213 
    2214        printf ("%5d ", (int)iter);
    2215        for (i = 0; i < (size_t)np; i++)
    2216          {
    2217            printf ("%10.3e ", gsl_vector_get (s->x, i));
    2218          }
    2219        printf ("f() = %7.3f size = %.3f\n", s->fval, size);
    2220      }
    2221    while (status == GSL_CONTINUE && iter < 100);
    2222 
    2223   for (i=0;i<(size_t)np;i++)
    2224     gsl_vector_set(par->x, i, gsl_vector_get(s->x, i));
    2225    //gsl_vector_free(par->x);
    2226    gsl_vector_free(ss);
    2227    gsl_multimin_fminimizer_free (s);
    2228716};
    2229717
    … …  
    2457945};
    2458946
    2459 /** Counts all cyclic bonds and returns their number.
    2460  * \note Hydrogen bonds can never by cyclic, thus no check for that
    2461  * \param *out output stream for debugging
    2462  * \return number opf cyclic bonds
    2463  */
    2464 int molecule::CountCyclicBonds(ofstream *out)
    2465 {
    2466   int No = 0;
    2467   int *MinimumRingSize = NULL;
    2468   MoleculeLeafClass *Subgraphs = NULL;
    2469   class StackClass<bond *> *BackEdgeStack = NULL;
    2470   bond *Binder = first;
    2471   if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
    2472     *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
    2473     Subgraphs = DepthFirstSearchAnalysis(out, BackEdgeStack);
    2474     while (Subgraphs->next != NULL) {
    2475       Subgraphs = Subgraphs->next;
    2476       delete(Subgraphs->previous);
    2477     }
    2478     delete(Subgraphs);
    2479     delete[](MinimumRingSize);
    2480   }
    2481   while(Binder->next != last) {
    2482     Binder = Binder->next;
    2483     if (Binder->Cyclic)
    2484       No++;
    2485   }
    2486   delete(BackEdgeStack);
    2487   return No;
    2488 };
    2489 /** Returns Shading as a char string.
    2490  * \param color the Shading
    2491  * \return string of the flag
    2492  */
    2493 string molecule::GetColor(enum Shading color)
    2494 {
    2495   switch(color) {
    2496     case white:
    2497       return "white";
    2498       break;
    2499     case lightgray:
    2500       return "lightgray";
    2501       break;
    2502     case darkgray:
    2503       return "darkgray";
    2504       break;
    2505     case black:
    2506       return "black";
    2507       break;
    2508     default:
    2509       return "uncolored";
    2510       break;
    2511   };
    2512 };
    2513947
    2514948
    … …  
    2539973};
    2540974
    2541 /** Creates an adjacency list of the molecule.
    2542  * We obtain an outside file with the indices of atoms which are bondmembers.
    2543  */
    2544 void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)
    2545 {
    2546 
    2547   // 1 We will parse bonds out of the dbond file created by tremolo.
    2548       int atom1, atom2, temp;
    2549       atom *Walker, *OtherWalker;
    2550 
    2551           if (!input)
    2552           {
    2553             cout << Verbose(1) << "Opening silica failed \n";
    2554           };
    2555 
    2556       *input >> ws >> atom1;
    2557       *input >> ws >> atom2;
    2558           cout << Verbose(1) << "Scanning file\n";
    2559           while (!input->eof()) // Check whether we read everything already
    2560           {
    2561         *input >> ws >> atom1;
    2562         *input >> ws >> atom2;
    2563             if(atom2<atom1) //Sort indices of atoms in order
    2564             {
    2565               temp=atom1;
    2566               atom1=atom2;
    2567               atom2=temp;
    2568             };
    2569 
    2570             Walker=start;
    2571             while(Walker-> nr != atom1) // Find atom corresponding to first index
    2572             {
    2573               Walker = Walker->next;
    2574             };
    2575             OtherWalker = Walker->next;
    2576             while(OtherWalker->nr != atom2) // Find atom corresponding to second index
    2577             {
    2578               OtherWalker= OtherWalker->next;
    2579             };
    2580             AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
    2581 
    2582           }
    2583 
    2584           CreateListOfBondsPerAtom(out);
    2585 
    2586 };
    2587 
    2588 
    2589 /** Creates an adjacency list of the molecule.
    2590  * Generally, we use the CSD approach to bond recognition, that is the the distance
    2591  * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
    2592  * a threshold t = 0.4 Angstroem.
    2593  * To make it O(N log N) the function uses the linked-cell technique as follows:
    2594  * The procedure is step-wise:
    2595  *  -# Remove every bond in list
    2596  *  -# Count the atoms in the molecule with CountAtoms()
    2597  *  -# partition cell into smaller linked cells of size \a bonddistance
    2598  *  -# put each atom into its corresponding cell
    2599  *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
    2600  *  -# create the list of bonds via CreateListOfBondsPerAtom()
    2601  *  -# correct the bond degree iteratively (single->double->triple bond)
    2602  *  -# finally print the bond list to \a *out if desired
    2603  * \param *out out stream for printing the matrix, NULL if no output
    2604  * \param bonddistance length of linked cells (i.e. maximum minimal length checked)
    2605  * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
    2606  */
    2607 void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
    2608 {
    2609 
    2610   atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
    2611   int No, NoBonds, CandidateBondNo;
    2612   int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
    2613   molecule **CellList;
    2614   double distance, MinDistance, MaxDistance;
    2615   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    2616   Vector x;
    2617   int FalseBondDegree = 0;
    2618 
    2619   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
    2620   *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
    2621   // remove every bond from the list
    2622   if ((first->next != last) && (last->previous != first)) {  // there are bonds present
    2623     cleanup(first,last);
    2624   }
    2625 
    2626   // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
    2627   CountAtoms(out);
    2628   *out << Verbose(1) << "AtomCount " << AtomCount << "." << endl;
    2629 
    2630   if (AtomCount != 0) {
    2631     // 1. find divisor for each axis, such that a sphere with radius of at least bonddistance can be placed into each cell
    2632     j=-1;
    2633     for (int i=0;i<NDIM;i++) {
    2634       j += i+1;
    2635       divisor[i] = (int)floor(cell_size[j]/bonddistance); // take smaller value such that size of linked cell is at least bonddistance
    2636       //*out << Verbose(1) << "divisor[" << i << "]  = " << divisor[i] << "." << endl;
    2637     }
    2638     // 2a. allocate memory for the cell list
    2639     NumberCells = divisor[0]*divisor[1]*divisor[2];
    2640     *out << Verbose(1) << "Allocating " << NumberCells << " cells." << endl;
    2641     CellList = Malloc<molecule*>(NumberCells, "molecule::CreateAdjacencyList - ** CellList");
    2642     for (int i=NumberCells;i--;)
    2643       CellList[i] = NULL;
    2644 
    2645     // 2b. put all atoms into its corresponding list
    2646     Walker = start;
    2647     while(Walker->next != end) {
    2648       Walker = Walker->next;
    2649       //*out << Verbose(1) << "Current atom is " << *Walker << " with coordinates ";
    2650       //Walker->x.Output(out);
    2651       //*out << "." << endl;
    2652       // compute the cell by the atom's coordinates
    2653       j=-1;
    2654       for (int i=0;i<NDIM;i++) {
    2655         j += i+1;
    2656         x.CopyVector(&(Walker->x));
    2657         x.KeepPeriodic(out, matrix);
    2658         n[i] = (int)floor(x.x[i]/cell_size[j]*(double)divisor[i]);
    2659       }
    2660       index = n[2] + (n[1] + n[0] * divisor[1]) * divisor[2];
    2661       //*out << Verbose(1) << "Atom " << *Walker << " goes into cell number [" << n[0] << "," << n[1] << "," << n[2] << "] = " << index << "." << endl;
    2662       // add copy atom to this cell
    2663       if (CellList[index] == NULL)  // allocate molecule if not done
    2664         CellList[index] = new molecule(elemente);
    2665       OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
    2666       //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
    2667     }
    2668     //for (int i=0;i<NumberCells;i++)
    2669       //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
    2670 
    2671 
    2672     // 3a. go through every cell
    2673     for (N[0]=divisor[0];N[0]--;)
    2674       for (N[1]=divisor[1];N[1]--;)
    2675         for (N[2]=divisor[2];N[2]--;) {
    2676           Index = N[2] + (N[1] + N[0] * divisor[1]) * divisor[2];
    2677           if (CellList[Index] != NULL) { // if there atoms in this cell
    2678             //*out << Verbose(1) << "Current cell is " << Index << "." << endl;
    2679             // 3b. for every atom therein
    2680             Walker = CellList[Index]->start;
    2681             while (Walker->next != CellList[Index]->end) {  // go through every atom
    2682               Walker = Walker->next;
    2683               //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
    2684               // 3c. check for possible bond between each atom in this and every one in the 27 cells
    2685               for (n[0]=-1;n[0]<=1;n[0]++)
    2686                 for (n[1]=-1;n[1]<=1;n[1]++)
    2687                   for (n[2]=-1;n[2]<=1;n[2]++) {
    2688                      // compute the index of this comparison cell and make it periodic
    2689                     index = ((N[2]+n[2]+divisor[2])%divisor[2]) + (((N[1]+n[1]+divisor[1])%divisor[1]) + ((N[0]+n[0]+divisor[0])%divisor[0]) * divisor[1]) * divisor[2];
    2690                     //*out << Verbose(1) << "Number of comparison cell is " << index << "." << endl;
    2691                     if (CellList[index] != NULL) {  // if there are any atoms in this cell
    2692                       OtherWalker = CellList[index]->start;
    2693                       while(OtherWalker->next != CellList[index]->end) {  // go through every atom in this cell
    2694                         OtherWalker = OtherWalker->next;
    2695                         //*out << Verbose(0) << "Current comparison atom is " << *OtherWalker << "." << endl;
    2696                         /// \todo periodic check is missing here!
    2697                         //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
    2698                         MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
    2699                         MinDistance *= (IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem;
    2700                         MaxDistance = MinDistance + BONDTHRESHOLD;
    2701                         MinDistance -= BONDTHRESHOLD;
    2702                         distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
    2703                         if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
    2704                           //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
    2705                           AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
    2706                         } else {
    2707                           //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
    2708                         }
    2709                       }
    2710                     }
    2711                   }
    2712             }
    2713           }
    2714         }
    2715 
    2716 
    2717 
    2718     // 4. free the cell again
    2719     for (int i=NumberCells;i--;)
    2720       if (CellList[i] != NULL) {
    2721         delete(CellList[i]);
    2722       }
    2723     Free(&CellList);
    2724 
    2725     // create the adjacency list per atom
    2726     CreateListOfBondsPerAtom(out);
    2727 
    2728     // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
    2729     // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
    2730     // preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
    2731     // double bonds as was expected.
    2732     if (BondCount != 0) {
    2733       NoCyclicBonds = 0;
    2734       *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
    2735       do {
    2736         No = 0; // No acts as breakup flag (if 1 we still continue)
    2737         Walker = start;
    2738         while (Walker->next != end) { // go through every atom
    2739           Walker = Walker->next;
    2740           // count valence of first partner
    2741           NoBonds = 0;
    2742           for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
    2743             NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
    2744           *out << Verbose(3) << "Walker " << *Walker << ": " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
    2745           if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
    2746             Candidate = NULL;
    2747             CandidateBondNo = -1;
    2748             for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
    2749               OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
    2750               // count valence of second partner
    2751               NoBonds = 0;
    2752               for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
    2753                 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
    2754               *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
    2755               if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
    2756                 if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
    2757                   Candidate = OtherWalker;
    2758                   CandidateBondNo = i;
    2759                   *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
    2760                 }
    2761               }
    2762             }
    2763             if ((Candidate != NULL) && (CandidateBondNo != -1)) {
    2764               ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
    2765               *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
    2766             } else
    2767               *out << Verbose(2) << "Could not find correct degree for atom " << *Walker << "." << endl;
    2768               FalseBondDegree++;
    2769           }
    2770         }
    2771       } while (No);
    2772     *out << " done." << endl;
    2773     } else
    2774       *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
    2775     *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
    2776 
    2777     // output bonds for debugging (if bond chain list was correctly installed)
    2778     *out << Verbose(1) << endl << "From contents of bond chain list:";
    2779     bond *Binder = first;
    2780     while(Binder->next != last) {
    2781       Binder = Binder->next;
    2782       *out << *Binder << "\t" << endl;
    2783     }
    2784     *out << endl;
    2785   } else
    2786     *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
    2787   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
    2788   Free(&matrix);
    2789 
    2790 };
    2791 
    2792 
    2793 
    2794 /** Performs a Depth-First search on this molecule.
    2795  * Marks bonds in molecule as cyclic, bridge, ... and atoms as
    2796  * articulations points, ...
    2797  * We use the algorithm from [Even, Graph Algorithms, p.62].
    2798  * \param *out output stream for debugging
    2799  * \param *&BackEdgeStack NULL pointer to StackClass with all the found back edges, allocated and filled on return
    2800  * \return list of each disconnected subgraph as an individual molecule class structure
    2801  */
    2802 MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack)
    2803 {
    2804   class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
    2805   BackEdgeStack = new StackClass<bond *> (BondCount);
    2806   MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
    2807   MoleculeLeafClass *LeafWalker = SubGraphs;
    2808   int CurrentGraphNr = 0, OldGraphNr;
    2809   int ComponentNumber = 0;
    2810   atom *Walker = NULL, *OtherAtom = NULL, *Root = start->next;
    2811   bond *Binder = NULL;
    2812   bool BackStepping = false;
    2813 
    2814   *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
    2815 
    2816   ResetAllBondsToUnused();
    2817   ResetAllAtomNumbers();
    2818   InitComponentNumbers();
    2819   BackEdgeStack->ClearStack();
    2820   while (Root != end) { // if there any atoms at all
    2821     // (1) mark all edges unused, empty stack, set atom->GraphNr = 0 for all
    2822     AtomStack->ClearStack();
    2823 
    2824     // put into new subgraph molecule and add this to list of subgraphs
    2825     LeafWalker = new MoleculeLeafClass(LeafWalker);
    2826     LeafWalker->Leaf = new molecule(elemente);
    2827     LeafWalker->Leaf->AddCopyAtom(Root);
    2828 
    2829     OldGraphNr = CurrentGraphNr;
    2830     Walker = Root;
    2831     do { // (10)
    2832       do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom
    2833         if (!BackStepping) { // if we don't just return from (8)
    2834           Walker->GraphNr = CurrentGraphNr;
    2835           Walker->LowpointNr = CurrentGraphNr;
    2836           *out << Verbose(1) << "Setting Walker[" << Walker->Name << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "." << endl;
    2837           AtomStack->Push(Walker);
    2838           CurrentGraphNr++;
    2839         }
    2840         do { // (3) if Walker has no unused egdes, go to (5)
    2841           BackStepping = false; // reset backstepping flag for (8)
    2842           if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
    2843             Binder = FindNextUnused(Walker);
    2844           if (Binder == NULL)
    2845             break;
    2846           *out << Verbose(2) << "Current Unused Bond is " << *Binder << "." << endl;
    2847           // (4) Mark Binder used, ...
    2848           Binder->MarkUsed(black);
    2849           OtherAtom = Binder->GetOtherAtom(Walker);
    2850           *out << Verbose(2) << "(4) OtherAtom is " << OtherAtom->Name << "." << endl;
    2851           if (OtherAtom->GraphNr != -1) {
    2852             // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
    2853             Binder->Type = BackEdge;
    2854             BackEdgeStack->Push(Binder);
    2855             Walker->LowpointNr = ( Walker->LowpointNr < OtherAtom->GraphNr ) ? Walker->LowpointNr : OtherAtom->GraphNr;
    2856             *out << Verbose(3) << "(4a) Visited: Setting Lowpoint of Walker[" << Walker->Name << "] to " << Walker->LowpointNr << "." << endl;
    2857           } else {
    2858             // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
    2859             Binder->Type = TreeEdge;
    2860             OtherAtom->Ancestor = Walker;
    2861             Walker = OtherAtom;
    2862             *out << Verbose(3) << "(4b) Not Visited: OtherAtom[" << OtherAtom->Name << "]'s Ancestor is now " << OtherAtom->Ancestor->Name << ", Walker is OtherAtom " << OtherAtom->Name << "." << endl;
    2863             break;
    2864           }
    2865           Binder = NULL;
    2866         } while (1);  // (3)
    2867         if (Binder == NULL) {
    2868           *out << Verbose(2) << "No more Unused Bonds." << endl;
    2869           break;
    2870         } else
    2871           Binder = NULL;
    2872       } while (1);  // (2)
    2873 
    2874       // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
    2875       if ((Walker == Root) && (Binder == NULL))
    2876         break;
    2877 
    2878       // (5) if Ancestor of Walker is ...
    2879       *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
    2880       if (Walker->Ancestor->GraphNr != Root->GraphNr) {
    2881         // (6)  (Ancestor of Walker is not Root)
    2882         if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
    2883           // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
    2884           Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
    2885           *out << Verbose(2) << "(6) Setting Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "." << endl;
    2886         } else {
    2887           // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
    2888           Walker->Ancestor->SeparationVertex = true;
    2889           *out << Verbose(2) << "(7) Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s is a separating vertex, creating component." << endl;
    2890           SetNextComponentNumber(Walker->Ancestor, ComponentNumber);
    2891           *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Ancestor's Compont is " << ComponentNumber << "." << endl;
    2892           SetNextComponentNumber(Walker, ComponentNumber);
    2893           *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Compont is " << ComponentNumber << "." << endl;
    2894           do {
    2895             OtherAtom = AtomStack->PopLast();
    2896             LeafWalker->Leaf->AddCopyAtom(OtherAtom);
    2897             SetNextComponentNumber(OtherAtom, ComponentNumber);
    2898             *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << ComponentNumber << "." << endl;
    2899           } while (OtherAtom != Walker);
    2900           ComponentNumber++;
    2901         }
    2902         // (8) Walker becomes its Ancestor, go to (3)
    2903         *out << Verbose(2) << "(8) Walker[" << Walker->Name << "] is now its Ancestor " << Walker->Ancestor->Name << ", backstepping. " << endl;
    2904         Walker = Walker->Ancestor;
    2905         BackStepping = true;
    2906       }
    2907       if (!BackStepping) {  // coming from (8) want to go to (3)
    2908         // (9) remove all from stack till Walker (including), these and Root form a component
    2909         AtomStack->Output(out);
    2910         SetNextComponentNumber(Root, ComponentNumber);
    2911         *out << Verbose(3) << "(9) Root[" << Root->Name << "]'s Component is " << ComponentNumber << "." << endl;
    2912         SetNextComponentNumber(Walker, ComponentNumber);
    2913         *out << Verbose(3) << "(9) Walker[" << Walker->Name << "]'s Component is " << ComponentNumber << "." << endl;
    2914         do {
    2915           OtherAtom = AtomStack->PopLast();
    2916           LeafWalker->Leaf->AddCopyAtom(OtherAtom);
    2917           SetNextComponentNumber(OtherAtom, ComponentNumber);
    2918           *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << ComponentNumber << "." << endl;
    2919         } while (OtherAtom != Walker);
    2920         ComponentNumber++;
    2921 
    2922         // (11) Root is separation vertex,  set Walker to Root and go to (4)
    2923         Walker = Root;
    2924         Binder = FindNextUnused(Walker);
    2925         *out << Verbose(1) << "(10) Walker is Root[" << Root->Name << "], next Unused Bond is " << Binder << "." << endl;
    2926         if (Binder != NULL) { // Root is separation vertex
    2927           *out << Verbose(1) << "(11) Root is a separation vertex." << endl;
    2928           Walker->SeparationVertex = true;
    2929         }
    2930       }
    2931     } while ((BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
    2932 
    2933     // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
    2934     *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
    2935     LeafWalker->Leaf->Output(out);
    2936     *out << endl;
    2937 
    2938     // step on to next root
    2939     while ((Root != end) && (Root->GraphNr != -1)) {
    2940       //*out << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
    2941       if (Root->GraphNr != -1) // if already discovered, step on
    2942         Root = Root->next;
    2943     }
    2944   }
    2945   // set cyclic bond criterium on "same LP" basis
    2946   Binder = first;
    2947   while(Binder->next != last) {
    2948     Binder = Binder->next;
    2949     if (Binder->rightatom->LowpointNr == Binder->leftatom->LowpointNr) { // cyclic ??
    2950       Binder->Cyclic = true;
    2951       NoCyclicBonds++;
    2952     }
    2953   }
    2954 
    2955 
    2956   *out << Verbose(1) << "Final graph info for each atom is:" << endl;
    2957   Walker = start;
    2958   while (Walker->next != end) {
    2959     Walker = Walker->next;
    2960     *out << Verbose(2) << "Atom " << Walker->Name << " is " << ((Walker->SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
    2961     OutputComponentNumber(out, Walker);
    2962     *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
    2963   }
    2964 
    2965   *out << Verbose(1) << "Final graph info for each bond is:" << endl;
    2966   Binder = first;
    2967   while(Binder->next != last) {
    2968     Binder = Binder->next;
    2969     *out << Verbose(2) << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
    2970     *out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
    2971     OutputComponentNumber(out, Binder->leftatom);
    2972     *out << " ===  ";
    2973     *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
    2974     OutputComponentNumber(out, Binder->rightatom);
    2975     *out << ">." << endl;
    2976     if (Binder->Cyclic) // cyclic ??
    2977       *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
    2978   }
    2979 
    2980   // free all and exit
    2981   delete(AtomStack);
    2982   *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
    2983   return SubGraphs;
    2984 };
    2985 
    2986 /** Analyses the cycles found and returns minimum of all cycle lengths.
    2987  * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root,
    2988  * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
    2989  * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
    2990  * as cyclic and print out the cycles.
    2991  * \param *out output stream for debugging
    2992  * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
    2993  * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
    2994  * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
    2995  */
    2996 void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *  BackEdgeStack, int *&MinimumRingSize)
    2997 {
    2998   atom **PredecessorList = Malloc<atom*>(AtomCount, "molecule::CyclicStructureAnalysis: **PredecessorList");
    2999   int *ShortestPathList = Malloc<int>(AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
    3000   enum Shading *ColorList = Malloc<enum Shading>(AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
    3001   class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring
    3002   class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
    3003   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
    3004   bond *Binder = NULL, *BackEdge = NULL;
    3005   int RingSize, NumCycles, MinRingSize = -1;
    3006 
    3007   // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
    3008   for (int i=AtomCount;i--;) {
    3009     PredecessorList[i] = NULL;
    3010     ShortestPathList[i] = -1;
    3011     ColorList[i] = white;
    3012   }
    3013 
    3014   *out << Verbose(1) << "Back edge list - ";
    3015   BackEdgeStack->Output(out);
    3016 
    3017   *out << Verbose(1) << "Analysing cycles ... " << endl;
    3018   NumCycles = 0;
    3019   while (!BackEdgeStack->IsEmpty()) {
    3020     BackEdge = BackEdgeStack->PopFirst();
    3021     // this is the target
    3022     Root = BackEdge->leftatom;
    3023     // this is the source point
    3024     Walker = BackEdge->rightatom;
    3025     ShortestPathList[Walker->nr] = 0;
    3026     BFSStack->ClearStack();  // start with empty BFS stack
    3027     BFSStack->Push(Walker);
    3028     TouchedStack->Push(Walker);
    3029     *out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
    3030     OtherAtom = NULL;
    3031     do {  // look for Root
    3032       Walker = BFSStack->PopFirst();
    3033       *out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
    3034       for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
    3035         Binder = ListOfBondsPerAtom[Walker->nr][i];
    3036         if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
    3037           OtherAtom = Binder->GetOtherAtom(Walker);
    3038 #ifdef ADDHYDROGEN
    3039           if (OtherAtom->type->Z != 1) {
    3040 #endif
    3041             *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
    3042             if (ColorList[OtherAtom->nr] == white) {
    3043               TouchedStack->Push(OtherAtom);
    3044               ColorList[OtherAtom->nr] = lightgray;
    3045               PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
    3046               ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
    3047               *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
    3048               //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
    3049                 *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
    3050                 BFSStack->Push(OtherAtom);
    3051               //}
    3052             } else {
    3053               *out << Verbose(3) << "Not Adding, has already been visited." << endl;
    3054             }
    3055             if (OtherAtom == Root)
    3056               break;
    3057 #ifdef ADDHYDROGEN
    3058           } else {
    3059             *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
    3060             ColorList[OtherAtom->nr] = black;
    3061           }
    3062 #endif
    3063         } else {
    3064           *out << Verbose(2) << "Bond " << *Binder << " not Visiting, is the back edge." << endl;
    3065         }
    3066       }
    3067       ColorList[Walker->nr] = black;
    3068       *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
    3069       if (OtherAtom == Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
    3070         // step through predecessor list
    3071         while (OtherAtom != BackEdge->rightatom) {
    3072           if (!OtherAtom->GetTrueFather()->IsCyclic)  // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
    3073             break;
    3074           else
    3075             OtherAtom = PredecessorList[OtherAtom->nr];
    3076         }
    3077         if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
    3078           *out << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl;\
    3079           do {
    3080             OtherAtom = TouchedStack->PopLast();
    3081             if (PredecessorList[OtherAtom->nr] == Walker) {
    3082               *out << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl;
    3083               PredecessorList[OtherAtom->nr] = NULL;
    3084               ShortestPathList[OtherAtom->nr] = -1;
    3085               ColorList[OtherAtom->nr] = white;
    3086               BFSStack->RemoveItem(OtherAtom);
    3087             }
    3088           } while ((!TouchedStack->IsEmpty()) && (PredecessorList[OtherAtom->nr] == NULL));
    3089           TouchedStack->Push(OtherAtom);  // last was wrongly popped
    3090           OtherAtom = BackEdge->rightatom; // set to not Root
    3091         } else
    3092           OtherAtom = Root;
    3093       }
    3094     } while ((!BFSStack->IsEmpty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
    3095 
    3096     if (OtherAtom == Root) {
    3097       // now climb back the predecessor list and thus find the cycle members
    3098       NumCycles++;
    3099       RingSize = 1;
    3100       Root->GetTrueFather()->IsCyclic = true;
    3101       *out << Verbose(1) << "Found ring contains: ";
    3102       Walker = Root;
    3103       while (Walker != BackEdge->rightatom) {
    3104         *out << Walker->Name << " <-> ";
    3105         Walker = PredecessorList[Walker->nr];
    3106         Walker->GetTrueFather()->IsCyclic = true;
    3107         RingSize++;
    3108       }
    3109       *out << Walker->Name << "  with a length of " << RingSize << "." << endl << endl;
    3110       // walk through all and set MinimumRingSize
    3111       Walker = Root;
    3112       MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
    3113       while (Walker != BackEdge->rightatom) {
    3114         Walker = PredecessorList[Walker->nr];
    3115         if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
    3116           MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
    3117       }
    3118       if ((RingSize < MinRingSize) || (MinRingSize == -1))
    3119         MinRingSize = RingSize;
    3120     } else {
    3121       *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
    3122     }
    3123 
    3124     // now clean the lists
    3125     while (!TouchedStack->IsEmpty()){
    3126       Walker = TouchedStack->PopFirst();
    3127       PredecessorList[Walker->nr] = NULL;
    3128       ShortestPathList[Walker->nr] = -1;
    3129       ColorList[Walker->nr] = white;
    3130     }
    3131   }
    3132   if (MinRingSize != -1) {
    3133     // go over all atoms
    3134     Root = start;
    3135     while(Root->next != end) {
    3136       Root = Root->next;
    3137 
    3138       if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
    3139         Walker = Root;
    3140         ShortestPathList[Walker->nr] = 0;
    3141         BFSStack->ClearStack();  // start with empty BFS stack
    3142         BFSStack->Push(Walker);
    3143         TouchedStack->Push(Walker);
    3144         //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
    3145         OtherAtom = Walker;
    3146         while (OtherAtom != NULL) {  // look for Root
    3147           Walker = BFSStack->PopFirst();
    3148           //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
    3149           for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
    3150             Binder = ListOfBondsPerAtom[Walker->nr][i];
    3151             if ((Binder != BackEdge) || (NumberOfBondsPerAtom[Walker->nr] == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
    3152               OtherAtom = Binder->GetOtherAtom(Walker);
    3153               //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
    3154               if (ColorList[OtherAtom->nr] == white) {
    3155                 TouchedStack->Push(OtherAtom);
    3156                 ColorList[OtherAtom->nr] = lightgray;
    3157                 PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
    3158                 ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
    3159                 //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
    3160                 if (OtherAtom->GetTrueFather()->IsCyclic) { // if the other atom is connected to a ring
    3161                   MinimumRingSize[Root->GetTrueFather()->nr] = ShortestPathList[OtherAtom->nr]+MinimumRingSize[OtherAtom->GetTrueFather()->nr];
    3162                   OtherAtom = NULL; //break;
    3163                   break;
    3164                 } else
    3165                   BFSStack->Push(OtherAtom);
    3166               } else {
    3167                 //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
    3168               }
    3169             } else {
    3170               //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
    3171             }
    3172           }
    3173           ColorList[Walker->nr] = black;
    3174           //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
    3175         }
    3176 
    3177         // now clean the lists
    3178         while (!TouchedStack->IsEmpty()){
    3179           Walker = TouchedStack->PopFirst();
    3180           PredecessorList[Walker->nr] = NULL;
    3181           ShortestPathList[Walker->nr] = -1;
    3182           ColorList[Walker->nr] = white;
    3183         }
    3184       }
    3185       *out << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl;
    3186     }
    3187     *out << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl;
    3188   } else
    3189     *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
    3190 
    3191   Free(&PredecessorList);
    3192   Free(&ShortestPathList);
    3193   Free(&ColorList);
    3194   delete(BFSStack);
    3195 };
    3196 
    3197 /** Sets the next component number.
    3198  * This is O(N) as the number of bonds per atom is bound.
    3199  * \param *vertex atom whose next atom::*ComponentNr is to be set
    3200  * \param nr number to use
    3201  */
    3202 void molecule::SetNextComponentNumber(atom *vertex, int nr)
    3203 {
    3204   int i=0;
    3205   if (vertex != NULL) {
    3206     for(;i<NumberOfBondsPerAtom[vertex->nr];i++) {
    3207       if (vertex->ComponentNr[i] == -1) {   // check if not yet used
    3208         vertex->ComponentNr[i] = nr;
    3209         break;
    3210       }
    3211       else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
    3212         break;  // breaking here will not cause error!
    3213     }
    3214     if (i == NumberOfBondsPerAtom[vertex->nr])
    3215       cerr << "Error: All Component entries are already occupied!" << endl;
    3216   } else
    3217       cerr << "Error: Given vertex is NULL!" << endl;
    3218 };
    3219 
    3220 /** Output a list of flags, stating whether the bond was visited or not.
    3221  * \param *out output stream for debugging
    3222  */
    3223 void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
    3224 {
    3225   for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
    3226     *out << vertex->ComponentNr[i] << "  ";
    3227 };
    3228 
    3229 /** Allocates memory for all atom::*ComponentNr in this molecule and sets each entry to -1.
    3230  */
    3231 void molecule::InitComponentNumbers()
    3232 {
    3233   atom *Walker = start;
    3234   while(Walker->next != end) {
    3235     Walker = Walker->next;
    3236     if (Walker->ComponentNr != NULL)
    3237       Free(&Walker->ComponentNr);
    3238     Walker->ComponentNr = Malloc<int>(NumberOfBondsPerAtom[Walker->nr], "molecule::InitComponentNumbers: *Walker->ComponentNr");
    3239     for (int i=NumberOfBondsPerAtom[Walker->nr];i--;)
    3240       Walker->ComponentNr[i] = -1;
    3241   }
    3242 };
    3243 
    3244 /** Returns next unused bond for this atom \a *vertex or NULL of none exists.
    3245  * \param *vertex atom to regard
    3246  * \return bond class or NULL
    3247  */
    3248 bond * molecule::FindNextUnused(atom *vertex)
    3249 {
    3250   for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
    3251     if (ListOfBondsPerAtom[vertex->nr][i]->IsUsed() == white)
    3252       return(ListOfBondsPerAtom[vertex->nr][i]);
    3253   return NULL;
    3254 };
    3255 
    3256 /** Resets bond::Used flag of all bonds in this molecule.
    3257  * \return true - success, false - -failure
    3258  */
    3259 void molecule::ResetAllBondsToUnused()
    3260 {
    3261   bond *Binder = first;
    3262   while (Binder->next != last) {
    3263     Binder = Binder->next;
    3264     Binder->ResetUsed();
    3265   }
    3266 };
    3267 
    3268 /** Resets atom::nr to -1 of all atoms in this molecule.
    3269  */
    3270 void molecule::ResetAllAtomNumbers()
    3271 {
    3272   atom *Walker = start;
    3273   while (Walker->next != end) {
    3274     Walker = Walker->next;
    3275     Walker->GraphNr  = -1;
    3276   }
    3277 };
    3278 
    3279 /** Output a list of flags, stating whether the bond was visited or not.
    3280  * \param *out output stream for debugging
    3281  * \param *list
    3282  */
    3283 void OutputAlreadyVisited(ofstream *out, int *list)
    3284 {
    3285   *out << Verbose(4) << "Already Visited Bonds:\t";
    3286   for(int i=1;i<=list[0];i++) *out << Verbose(0) << list[i] << "  ";
    3287   *out << endl;
    3288 };
    3289 
    3290 /** Estimates by educated guessing (using upper limit) the expected number of fragments.
    3291  * The upper limit is
    3292  * \f[
    3293  *  n = N \cdot C^k
    3294  * \f]
    3295  * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
    3296  * \param *out output stream for debugging
    3297  * \param order bond order k
    3298  * \return number n of fragments
    3299  */
    3300 int molecule::GuesstimateFragmentCount(ofstream *out, int order)
    3301 {
    3302   int c = 0;
    3303   int FragmentCount;
    3304   // get maximum bond degree
    3305   atom *Walker = start;
    3306   while (Walker->next != end) {
    3307     Walker = Walker->next;
    3308     c = (NumberOfBondsPerAtom[Walker->nr] > c) ? NumberOfBondsPerAtom[Walker->nr] : c;
    3309   }
    3310   FragmentCount = NoNonHydrogen*(1 << (c*order));
    3311   *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
    3312   return FragmentCount;
    3313 };
    3314 
    3315 /** Scans a single line for number and puts them into \a KeySet.
    3316  * \param *out output stream for debugging
    3317  * \param *buffer buffer to scan
    3318  * \param &CurrentSet filled KeySet on return
    3319  * \return true - at least one valid atom id parsed, false - CurrentSet is empty
    3320  */
    3321 bool molecule::ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
    3322 {
    3323   stringstream line;
    3324   int AtomNr;
    3325   int status = 0;
    3326 
    3327   line.str(buffer);
    3328   while (!line.eof()) {
    3329     line >> AtomNr;
    3330     if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
    3331       CurrentSet.insert(AtomNr);  // insert at end, hence in same order as in file!
    3332       status++;
    3333     } // else it's "-1" or else and thus must not be added
    3334   }
    3335   *out << Verbose(1) << "The scanned KeySet is ";
    3336   for(KeySet::iterator runner = CurrentSet.begin(); runner != CurrentSet.end(); runner++) {
    3337     *out << (*runner) << "\t";
    3338   }
    3339   *out << endl;
    3340   return (status != 0);
    3341 };
    3342 
    3343 /** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
    3344  * Does two-pass scanning:
    3345  * -# Scans the keyset file and initialises a temporary graph
    3346  * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
    3347  * Finally, the temporary graph is inserted into the given \a FragmentList for return.
    3348  * \param *out output stream for debugging
    3349  * \param *path path to file
    3350  * \param *FragmentList empty, filled on return
    3351  * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
    3352  */
    3353 bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
    3354 {
    3355   bool status = true;
    3356   ifstream InputFile;
    3357   stringstream line;
    3358   GraphTestPair testGraphInsert;
    3359   int NumberOfFragments = 0;
    3360   double TEFactor;
    3361   char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
    3362 
    3363   if (FragmentList == NULL) { // check list pointer
    3364     FragmentList = new Graph;
    3365   }
    3366 
    3367   // 1st pass: open file and read
    3368   *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
    3369   sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
    3370   InputFile.open(filename);
    3371   if (InputFile != NULL) {
    3372     // each line represents a new fragment
    3373     char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - *buffer");
    3374     // 1. parse keysets and insert into temp. graph
    3375     while (!InputFile.eof()) {
    3376       InputFile.getline(buffer, MAXSTRINGSIZE);
    3377       KeySet CurrentSet;
    3378       if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet))) {  // if at least one valid atom was added, write config
    3379         testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1)));  // store fragment number and current factor
    3380         if (!testGraphInsert.second) {
    3381           cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
    3382         }
    3383       }
    3384     }
    3385     // 2. Free and done
    3386     InputFile.close();
    3387     InputFile.clear();
    3388     Free(&buffer);
    3389     *out << Verbose(1) << "done." << endl;
    3390   } else {
    3391     *out << Verbose(1) << "File " << filename << " not found." << endl;
    3392     status = false;
    3393   }
    3394 
    3395   // 2nd pass: open TEFactors file and read
    3396   *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
    3397   sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, TEFACTORSFILE);
    3398   InputFile.open(filename);
    3399   if (InputFile != NULL) {
    3400     // 3. add found TEFactors to each keyset
    3401     NumberOfFragments = 0;
    3402     for(Graph::iterator runner = FragmentList->begin();runner != FragmentList->end(); runner++) {
    3403       if (!InputFile.eof()) {
    3404         InputFile >> TEFactor;
    3405         (*runner).second.second = TEFactor;
    3406         *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
    3407       } else {
    3408         status = false;
    3409         break;
    3410       }
    3411     }
    3412     // 4. Free and done
    3413     InputFile.close();
    3414     *out << Verbose(1) << "done." << endl;
    3415   } else {
    3416     *out << Verbose(1) << "File " << filename << " not found." << endl;
    3417     status = false;
    3418   }
    3419 
    3420   // free memory
    3421   Free(&filename);
    3422 
    3423   return status;
    3424 };
    3425 
    3426 /** Stores keysets and TEFactors to file.
    3427  * \param *out output stream for debugging
    3428  * \param KeySetList Graph with Keysets and factors
    3429  * \param *path path to file
    3430  * \return true - file written successfully, false - writing failed
    3431  */
    3432 bool molecule::StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
    3433 {
    3434   ofstream output;
    3435   bool status =  true;
    3436   string line;
    3437 
    3438   // open KeySet file
    3439   line = path;
    3440   line.append("/");
    3441   line += FRAGMENTPREFIX;
    3442   line += KEYSETFILE;
    3443   output.open(line.c_str(), ios::out);
    3444   *out << Verbose(1) << "Saving key sets of the total graph ... ";
    3445   if(output != NULL) {
    3446     for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
    3447       for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    3448         if (sprinter != (*runner).first.begin())
    3449           output << "\t";
    3450         output << *sprinter;
    3451       }
    3452       output << endl;
    3453     }
    3454     *out << "done." << endl;
    3455   } else {
    3456     cerr << "Unable to open " << line << " for writing keysets!" << endl;
    3457     status = false;
    3458   }
    3459   output.close();
    3460   output.clear();
    3461 
    3462   // open TEFactors file
    3463   line = path;
    3464   line.append("/");
    3465   line += FRAGMENTPREFIX;
    3466   line += TEFACTORSFILE;
    3467   output.open(line.c_str(), ios::out);
    3468   *out << Verbose(1) << "Saving TEFactors of the total graph ... ";
    3469   if(output != NULL) {
    3470     for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++)
    3471       output << (*runner).second.second << endl;
    3472     *out << Verbose(1) << "done." << endl;
    3473   } else {
    3474     *out << Verbose(1) << "failed to open " << line << "." << endl;
    3475     status = false;
    3476   }
    3477   output.close();
    3478 
    3479   return status;
    3480 };
    3481 
    3482 /** Storing the bond structure of a molecule to file.
    3483  * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
    3484  * \param *out output stream for debugging
    3485  * \param *path path to file
    3486  * \return true - file written successfully, false - writing failed
    3487  */
    3488 bool molecule::StoreAdjacencyToFile(ofstream *out, char *path)
    3489 {
    3490   ofstream AdjacencyFile;
    3491   atom *Walker = NULL;
    3492   stringstream line;
    3493   bool status = true;
    3494 
    3495   line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
    3496   AdjacencyFile.open(line.str().c_str(), ios::out);
    3497   *out << Verbose(1) << "Saving adjacency list ... ";
    3498   if (AdjacencyFile != NULL) {
    3499     Walker = start;
    3500     while(Walker->next != end) {
    3501       Walker = Walker->next;
    3502       AdjacencyFile << Walker->nr << "\t";
    3503       for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
    3504         AdjacencyFile << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
    3505       AdjacencyFile << endl;
    3506     }
    3507     AdjacencyFile.close();
    3508     *out << Verbose(1) << "done." << endl;
    3509   } else {
    3510     *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    3511     status = false;
    3512   }
    3513 
    3514   return status;
    3515 };
    3516 
    3517 /** Checks contents of adjacency file against bond structure in structure molecule.
    3518  * \param *out output stream for debugging
    3519  * \param *path path to file
    3520  * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
    3521  * \return true - structure is equal, false - not equivalence
    3522  */
    3523 bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
    3524 {
    3525   ifstream File;
    3526   stringstream filename;
    3527   bool status = true;
    3528   char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
    3529 
    3530   filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
    3531   File.open(filename.str().c_str(), ios::out);
    3532   *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
    3533   if (File != NULL) {
    3534     // allocate storage structure
    3535     int NonMatchNumber = 0;   // will number of atoms with differing bond structure
    3536     int *CurrentBonds = Malloc<int>(8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
    3537     int CurrentBondsOfAtom;
    3538 
    3539     // Parse the file line by line and count the bonds
    3540     while (!File.eof()) {
    3541       File.getline(buffer, MAXSTRINGSIZE);
    3542       stringstream line;
    3543       line.str(buffer);
    3544       int AtomNr = -1;
    3545       line >> AtomNr;
    3546       CurrentBondsOfAtom = -1; // we count one too far due to line end
    3547       // parse into structure
    3548       if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
    3549         while (!line.eof())
    3550           line >> CurrentBonds[ ++CurrentBondsOfAtom ];
    3551         // compare against present bonds
    3552         //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
    3553         if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
    3554           for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
    3555             int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
    3556             int j = 0;
    3557             for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
    3558             if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
    3559               ListOfAtoms[AtomNr] = NULL;
    3560               NonMatchNumber++;
    3561               status = false;
    3562               //out << "[" << id << "]\t";
    3563             } else {
    3564               //out << id << "\t";
    3565             }
    3566           }
    3567           //out << endl;
    3568         } else {
    3569           *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
    3570           status = false;
    3571         }
    3572       }
    3573     }
    3574     File.close();
    3575     File.clear();
    3576     if (status) { // if equal we parse the KeySetFile
    3577       *out << Verbose(1) << "done: Equal." << endl;
    3578       status = true;
    3579     } else
    3580       *out << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl;
    3581     Free(&CurrentBonds);
    3582   } else {
    3583     *out << Verbose(1) << "Adjacency file not found." << endl;
    3584     status = false;
    3585   }
    3586   *out << endl;
    3587   Free(&buffer);
    3588 
    3589   return status;
    3590 };
    3591 
    3592 /** Checks whether the OrderAtSite is still below \a Order at some site.
    3593  * \param *out output stream for debugging
    3594  * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
    3595  * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
    3596  * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
    3597  * \param *MinimumRingSize array of max. possible order to avoid loops
    3598  * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
    3599  * \return true - needs further fragmentation, false - does not need fragmentation
    3600  */
    3601 bool molecule::CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
    3602 {
    3603   atom *Walker = start;
    3604   bool status = false;
    3605   ifstream InputFile;
    3606 
    3607   // initialize mask list
    3608   for(int i=AtomCount;i--;)
    3609     AtomMask[i] = false;
    3610 
    3611   if (Order < 0) { // adaptive increase of BondOrder per site
    3612     if (AtomMask[AtomCount] == true)  // break after one step
    3613       return false;
    3614     // parse the EnergyPerFragment file
    3615     char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
    3616     sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
    3617     InputFile.open(buffer, ios::in);
    3618     if ((InputFile != NULL) && (GlobalKeySetList != NULL)) {
    3619       // transmorph graph keyset list into indexed KeySetList
    3620       map<int,KeySet> IndexKeySetList;
    3621       for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
    3622         IndexKeySetList.insert( pair<int,KeySet>(runner->second.first,runner->first) );
    3623       }
    3624       int lines = 0;
    3625       // count the number of lines, i.e. the number of fragments
    3626       InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
    3627       InputFile.getline(buffer, MAXSTRINGSIZE);
    3628       while(!InputFile.eof()) {
    3629         InputFile.getline(buffer, MAXSTRINGSIZE);
    3630         lines++;
    3631       }
    3632       //*out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
    3633       InputFile.clear();
    3634       InputFile.seekg(ios::beg);
    3635       map<int, pair<double,int> > AdaptiveCriteriaList;  // (Root No., (Value, Order)) !
    3636       int No, FragOrder;
    3637       double Value;
    3638       // each line represents a fragment root (Atom::nr) id and its energy contribution
    3639       InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
    3640       InputFile.getline(buffer, MAXSTRINGSIZE);
    3641       while(!InputFile.eof()) {
    3642         InputFile.getline(buffer, MAXSTRINGSIZE);
    3643         if (strlen(buffer) > 2) {
    3644           //*out << Verbose(2) << "Scanning: " << buffer << endl;
    3645           stringstream line(buffer);
    3646           line >> FragOrder;
    3647           line >> ws >> No;
    3648           line >> ws >> Value; // skip time entry
    3649           line >> ws >> Value;
    3650           No -= 1;  // indices start at 1 in file, not 0
    3651           //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
    3652 
    3653           // clean the list of those entries that have been superceded by higher order terms already
    3654           map<int,KeySet>::iterator marker = IndexKeySetList.find(No);    // find keyset to Frag No.
    3655           if (marker != IndexKeySetList.end()) {  // if found
    3656             Value *= 1 + MYEPSILON*(*((*marker).second.begin()));     // in case of equal energies this makes em not equal without changing anything actually
    3657             // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
    3658             pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
    3659             map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
    3660             if (!InsertedElement.second) { // this root is already present
    3661               if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
    3662                 //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
    3663                 {  // if value is smaller, update value and order
    3664                 (*PresentItem).second.first = fabs(Value);
    3665                 (*PresentItem).second.second = FragOrder;
    3666                 *out << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
    3667               } else {
    3668                 *out << Verbose(2) << "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl;
    3669               }
    3670             } else {
    3671               *out << Verbose(2) << "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
    3672             }
    3673           } else {
    3674             *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
    3675           }
    3676         }
    3677       }
    3678       // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
    3679       map<double, pair<int,int> > FinalRootCandidates;
    3680       *out << Verbose(1) << "Root candidate list is: " << endl;
    3681       for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList.begin(); runner != AdaptiveCriteriaList.end(); runner++) {
    3682         Walker = FindAtom((*runner).first);
    3683         if (Walker != NULL) {
    3684           //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
    3685           if (!Walker->MaxOrder) {
    3686             *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
    3687             FinalRootCandidates.insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
    3688           } else {
    3689             *out << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl;
    3690           }
    3691         } else {
    3692           cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
    3693         }
    3694       }
    3695       // pick the ones still below threshold and mark as to be adaptively updated
    3696       for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
    3697         No = (*runner).second.first;
    3698         Walker = FindAtom(No);
    3699         //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
    3700           *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
    3701           AtomMask[No] = true;
    3702           status = true;
    3703         //} else
    3704           //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
    3705       }
    3706       // close and done
    3707       InputFile.close();
    3708       InputFile.clear();
    3709     } else {
    3710       cerr << "Unable to parse " << buffer << " file, incrementing all." << endl;
    3711       while (Walker->next != end) {
    3712         Walker = Walker->next;
    3713     #ifdef ADDHYDROGEN
    3714         if (Walker->type->Z != 1) // skip hydrogen
    3715     #endif
    3716         {
    3717           AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
    3718           status = true;
    3719         }
    3720       }
    3721     }
    3722     Free(&buffer);
    3723     // pick a given number of highest values and set AtomMask
    3724   } else { // global increase of Bond Order
    3725     while (Walker->next != end) {
    3726       Walker = Walker->next;
    3727   #ifdef ADDHYDROGEN
    3728       if (Walker->type->Z != 1) // skip hydrogen
    3729   #endif
    3730       {
    3731         AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
    3732         if ((Order != 0) && (Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
    3733           status = true;
    3734       }
    3735     }
    3736     if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
    3737       status = true;
    3738 
    3739     if (!status) {
    3740       if (Order == 0)
    3741         *out << Verbose(1) << "Single stepping done." << endl;
    3742       else
    3743         *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
    3744     }
    3745   }
    3746 
    3747   // print atom mask for debugging
    3748   *out << "              ";
    3749   for(int i=0;i<AtomCount;i++)
    3750     *out << (i % 10);
    3751   *out << endl << "Atom mask is: ";
    3752   for(int i=0;i<AtomCount;i++)
    3753     *out << (AtomMask[i] ? "t" : "f");
    3754   *out << endl;
    3755 
    3756   return status;
    3757 };
    3758 
    3759 /** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
    3760  * \param *out output stream for debugging
    3761  * \param *&SortIndex Mapping array of size molecule::AtomCount
    3762  * \return true - success, false - failure of SortIndex alloc
    3763  * \todo do we really need this still as the IonType may appear in any order due to recent changes
    3764  */
    3765 bool molecule::CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex)
    3766 {
    3767   element *runner = elemente->start;
    3768   int AtomNo = 0;
    3769   atom *Walker = NULL;
    3770 
    3771   if (SortIndex != NULL) {
    3772     *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
    3773     return false;
    3774   }
    3775   SortIndex = Malloc<int>(AtomCount, "molecule::FragmentMolecule: *SortIndex");
    3776   for(int i=AtomCount;i--;)
    3777     SortIndex[i] = -1;
    3778   while (runner->next != elemente->end) { // go through every element
    3779     runner = runner->next;
    3780     if (ElementsInMolecule[runner->Z]) { // if this element got atoms
    3781       Walker = start;
    3782       while (Walker->next != end) { // go through every atom of this element
    3783         Walker = Walker->next;
    3784         if (Walker->type->Z == runner->Z) // if this atom fits to element
    3785           SortIndex[Walker->nr] = AtomNo++;
    3786       }
    3787     }
    3788   }
    3789   return true;
    3790 };
    3791 
    3792 /** Performs a many-body bond order analysis for a given bond order.
    3793  * -# parses adjacency, keysets and orderatsite files
    3794  * -# performs DFS to find connected subgraphs (to leave this in was a design decision: might be useful later)
    3795  * -# RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energ
    3796 y contribution", and that's why this consciously not done in the following loop)
    3797  * -# in a loop over all subgraphs
    3798  *  -# calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
    3799  *  -# creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)
    3800  * -# combines the generated molecule lists from all subgraphs
    3801  * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
    3802  * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
    3803  * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
    3804  * subgraph in the MoleculeListClass.
    3805  * \param *out output stream for debugging
    3806  * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
    3807  * \param *configuration configuration for writing config files for each fragment
    3808  * \return 1 - continue, 2 - stop (no fragmentation occured)
    3809  */
    3810 int molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
    3811 {
    3812   MoleculeListClass *BondFragments = NULL;
    3813   int *SortIndex = NULL;
    3814   int *MinimumRingSize = new int[AtomCount];
    3815   int FragmentCounter;
    3816   MoleculeLeafClass *MolecularWalker = NULL;
    3817   MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
    3818   fstream File;
    3819   bool FragmentationToDo = true;
    3820   class StackClass<bond *> *BackEdgeStack = NULL, *LocalBackEdgeStack = NULL;
    3821   bool CheckOrder = false;
    3822   Graph **FragmentList = NULL;
    3823   Graph *ParsedFragmentList = NULL;
    3824   Graph TotalGraph;     // graph with all keysets however local numbers
    3825   int TotalNumberOfKeySets = 0;
    3826   atom **ListOfAtoms = NULL;
    3827   atom ***ListOfLocalAtoms = NULL;
    3828   bool *AtomMask = NULL;
    3829 
    3830   *out << endl;
    3831 #ifdef ADDHYDROGEN
    3832   *out << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
    3833 #else
    3834   *out << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl;
    3835 #endif
    3836 
    3837   // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
    3838 
    3839   // ===== 1. Check whether bond structure is same as stored in files ====
    3840 
    3841   // fill the adjacency list
    3842   CreateListOfBondsPerAtom(out);
    3843 
    3844   // create lookup table for Atom::nr
    3845   FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);
    3846 
    3847   // === compare it with adjacency file ===
    3848   FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
    3849   Free(&ListOfAtoms);
    3850 
    3851   // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
    3852   Subgraphs = DepthFirstSearchAnalysis(out, BackEdgeStack);
    3853   // fill the bond structure of the individually stored subgraphs
    3854   Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
    3855   // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
    3856   for(int i=AtomCount;i--;)
    3857     MinimumRingSize[i] = AtomCount;
    3858   MolecularWalker = Subgraphs;
    3859   FragmentCounter = 0;
    3860   while (MolecularWalker->next != NULL) {
    3861     MolecularWalker = MolecularWalker->next;
    3862     *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    3863     LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
    3864 //    // check the list of local atoms for debugging
    3865 //    *out << Verbose(0) << "ListOfLocalAtoms for this subgraph is:" << endl;
    3866 //    for (int i=0;i<AtomCount;i++)
    3867 //      if (ListOfLocalAtoms[FragmentCounter][i] == NULL)
    3868 //        *out << "\tNULL";
    3869 //      else
    3870 //        *out << "\t" << ListOfLocalAtoms[FragmentCounter][i]->Name;
    3871     *out << Verbose(0) << "Gathering local back edges for subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    3872     MolecularWalker->Leaf->PickLocalBackEdges(out, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
    3873     *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    3874     MolecularWalker->Leaf->CyclicStructureAnalysis(out, LocalBackEdgeStack, MinimumRingSize);
    3875     *out << Verbose(0) << "Done with Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    3876     delete(LocalBackEdgeStack);
    3877   }
    3878  
    3879   // ===== 3. if structure still valid, parse key set file and others =====
    3880   FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
    3881 
    3882   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
    3883   FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
    3884  
    3885   // =================================== Begin of FRAGMENTATION ===============================
    3886   // ===== 6a. assign each keyset to its respective subgraph =====
    3887   Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
    3888 
    3889   // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
    3890   KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
    3891   AtomMask = new bool[AtomCount+1];
    3892   AtomMask[AtomCount] = false;
    3893   FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
    3894   while ((CheckOrder = CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
    3895     FragmentationToDo = FragmentationToDo || CheckOrder;
    3896     AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    3897     // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    3898     Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
    3899 
    3900     // ===== 7. fill the bond fragment list =====
    3901     FragmentCounter = 0;
    3902     MolecularWalker = Subgraphs;
    3903     while (MolecularWalker->next != NULL) {
    3904       MolecularWalker = MolecularWalker->next;
    3905       *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
    3906       //MolecularWalker->Leaf->OutputListOfBonds(out);  // output ListOfBondsPerAtom for debugging
    3907       if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
    3908         // call BOSSANOVA method
    3909         *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
    3910         MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
    3911       } else {
    3912         cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
    3913       }
    3914       FragmentCounter++;  // next fragment list
    3915     }
    3916   }
    3917   delete[](RootStack);
    3918   delete[](AtomMask);
    3919   delete(ParsedFragmentList);
    3920   delete[](MinimumRingSize);
    3921  
    3922 
    3923   // ==================================== End of FRAGMENTATION ============================================
    3924 
    3925   // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
    3926   Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
    3927 
    3928   // free subgraph memory again
    3929   FragmentCounter = 0;
    3930   if (Subgraphs != NULL) {
    3931     while (Subgraphs->next != NULL) {
    3932       Subgraphs = Subgraphs->next;
    3933       delete(FragmentList[FragmentCounter++]);
    3934       delete(Subgraphs->previous);
    3935     }
    3936     delete(Subgraphs);
    3937   }
    3938   Free(&FragmentList);
    3939 
    3940   // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
    3941   //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
    3942     // allocate memory for the pointer array and transmorph graphs into full molecular fragments
    3943     BondFragments = new MoleculeListClass();
    3944     int k=0;
    3945     for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
    3946       KeySet test = (*runner).first;
    3947       *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
    3948       BondFragments->insert(StoreFragmentFromKeySet(out, test, configuration));
    3949       k++;
    3950     }
    3951     *out << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl;
    3952 
    3953     // ===== 9. Save fragments' configuration and keyset files et al to disk ===
    3954     if (BondFragments->ListOfMolecules.size() != 0) {
    3955       // create the SortIndex from BFS labels to order in the config file
    3956       CreateMappingLabelsToConfigSequence(out, SortIndex);
    3957 
    3958       *out << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl;
    3959       if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
    3960         *out << Verbose(1) << "All configs written." << endl;
    3961       else
    3962         *out << Verbose(1) << "Some config writing failed." << endl;
    3963 
    3964       // store force index reference file
    3965       BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
    3966 
    3967       // store keysets file
    3968       StoreKeySetFile(out, TotalGraph, configuration->configpath);
    3969 
    3970       // store Adjacency file
    3971       StoreAdjacencyToFile(out, configuration->configpath);
    3972 
    3973       // store Hydrogen saturation correction file
    3974       BondFragments->AddHydrogenCorrection(out, configuration->configpath);
    3975 
    3976       // store adaptive orders into file
    3977       StoreOrderAtSiteFile(out, configuration->configpath);
    3978 
    3979       // restore orbital and Stop values
    3980       CalculateOrbitals(*configuration);
    3981 
    3982       // free memory for bond part
    3983       *out << Verbose(1) << "Freeing bond memory" << endl;
    3984       delete(FragmentList); // remove bond molecule from memory
    3985       Free(&SortIndex);
    3986     } else
    3987       *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
    3988   //} else
    3989   //  *out << Verbose(1) << "No fragments to store." << endl;
    3990   *out << Verbose(0) << "End of bond fragmentation." << endl;
    3991 
    3992   return ((int)(!FragmentationToDo)+1);    // 1 - continue, 2 - stop (no fragmentation occured)
    3993 };
    3994 
    3995 
    3996 /** Picks from a global stack with all back edges the ones in the fragment.
    3997  * \param *out output stream for debugging
    3998  * \param **ListOfLocalAtoms array of father atom::nr to local atom::nr (reverse of atom::father)
    3999  * \param *ReferenceStack stack with all the back egdes
    4000  * \param *LocalStack stack to be filled
    4001  * \return true - everything ok, false - ReferenceStack was empty
    4002  */
    4003 bool molecule::PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack)
    4004 {
    4005   bool status = true;
    4006   if (ReferenceStack->IsEmpty()) {
    4007     cerr << "ReferenceStack is empty!" << endl;
    4008     return false;
    4009   }
    4010   bond *Binder = ReferenceStack->PopFirst();
    4011   bond *FirstBond = Binder;   // mark the first bond, so that we don't loop through the stack indefinitely
    4012   atom *Walker = NULL, *OtherAtom = NULL;
    4013   ReferenceStack->Push(Binder);
    4014  
    4015   do {  // go through all bonds and push local ones
    4016     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
    4017     if (Walker != NULL) // if this Walker exists in the subgraph ...
    4018       for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
    4019         OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
    4020         if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
    4021           LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
    4022           *out << Verbose(3) << "Found local edge " << *(ListOfBondsPerAtom[Walker->nr][i]) << "." << endl;
    4023           break;
    4024         }
    4025       }
    4026     Binder = ReferenceStack->PopFirst();  // loop the stack for next item
    4027     *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
    4028     ReferenceStack->Push(Binder);
    4029   } while (FirstBond != Binder);
    4030 
    4031   return status;
    4032 };
    4033 
    4034 /** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
    4035  * Atoms not present in the file get "-1".
    4036  * \param *out output stream for debugging
    4037  * \param *path path to file ORDERATSITEFILE
    4038  * \return true - file writable, false - not writable
    4039  */
    4040 bool molecule::StoreOrderAtSiteFile(ofstream *out, char *path)
    4041 {
    4042   stringstream line;
    4043   ofstream file;
    4044 
    4045   line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
    4046   file.open(line.str().c_str());
    4047   *out << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
    4048   if (file != NULL) {
    4049     atom *Walker = start;
    4050     while (Walker->next != end) {
    4051       Walker = Walker->next;
    4052       file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << endl;
    4053       *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << "." << endl;
    4054     }
    4055     file.close();
    4056     *out << Verbose(1) << "done." << endl;
    4057     return true;
    4058   } else {
    4059     *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    4060     return false;
    4061   }
    4062 };
    4063 
    4064 /** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
    4065  * Atoms not present in the file get "0".
    4066  * \param *out output stream for debugging
    4067  * \param *path path to file ORDERATSITEFILEe
    4068  * \return true - file found and scanned, false - file not found
    4069  * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
    4070  */
    4071 bool molecule::ParseOrderAtSiteFromFile(ofstream *out, char *path)
    4072 {
    4073   unsigned char *OrderArray = Malloc<unsigned char>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
    4074   bool *MaxArray = Malloc<bool>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
    4075   bool status;
    4076   int AtomNr, value;
    4077   stringstream line;
    4078   ifstream file;
    4079 
    4080   *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
    4081   for(int i=AtomCount;i--;)
    4082     OrderArray[i] = 0;
    4083   line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
    4084   file.open(line.str().c_str());
    4085   if (file != NULL) {
    4086     for (int i=AtomCount;i--;) { // initialise with 0
    4087       OrderArray[i] = 0;
    4088       MaxArray[i] = 0;
    4089     }
    4090     while (!file.eof()) { // parse from file
    4091       AtomNr = -1;
    4092       file >> AtomNr;
    4093       if (AtomNr != -1) {   // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
    4094         file >> value;
    4095         OrderArray[AtomNr] = value;
    4096         file >> value;
    4097         MaxArray[AtomNr] = value;
    4098         //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << "." << endl;
    4099       }
    4100     }
    4101     atom *Walker = start;
    4102     while (Walker->next != end) { // fill into atom classes
    4103       Walker = Walker->next;
    4104       Walker->AdaptiveOrder = OrderArray[Walker->nr];
    4105       Walker->MaxOrder = MaxArray[Walker->nr];
    4106       *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl;
    4107     }
    4108     file.close();
    4109     *out << Verbose(1) << "done." << endl;
    4110     status = true;
    4111   } else {
    4112     *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    4113     status = false;
    4114   }
    4115   Free(&OrderArray);
    4116   Free(&MaxArray);
    4117 
    4118   *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
    4119   return status;
    4120 };
    4121975
    4122976/** Creates an 2d array of pointer with an entry for each atom and each bond it has.
    … …  
    41851039};
    41861040
    4187 /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
    4188  * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
    4189  * white and putting into queue.
    4190  * \param *out output stream for debugging
    4191  * \param *Mol Molecule class to add atoms to
    4192  * \param **AddedAtomList list with added atom pointers, index is atom father's number
    4193  * \param **AddedBondList list with added bond pointers, index is bond father's number
    4194  * \param *Root root vertex for BFS
    4195  * \param *Bond bond not to look beyond
    4196  * \param BondOrder maximum distance for vertices to add
    4197  * \param IsAngstroem lengths are in angstroem or bohrradii
    4198  */
    4199 void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
    4200 {
    4201   atom **PredecessorList = Malloc<atom*>(AtomCount, "molecule::BreadthFirstSearchAdd: **PredecessorList");
    4202   int *ShortestPathList = Malloc<int>(AtomCount, "molecule::BreadthFirstSearchAdd: *ShortestPathList");
    4203   enum Shading *ColorList = Malloc<enum Shading>(AtomCount, "molecule::BreadthFirstSearchAdd: *ColorList");
    4204   class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
    4205   atom *Walker = NULL, *OtherAtom = NULL;
    4206   bond *Binder = NULL;
    4207 
    4208   // add Root if not done yet
    4209   AtomStack->ClearStack();
    4210   if (AddedAtomList[Root->nr] == NULL)  // add Root if not yet present
    4211     AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
    4212   AtomStack->Push(Root);
    4213 
    4214   // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
    4215   for (int i=AtomCount;i--;) {
    4216     PredecessorList[i] = NULL;
    4217     ShortestPathList[i] = -1;
    4218     if (AddedAtomList[i] != NULL) // mark already present atoms (i.e. Root and maybe others) as visited
    4219       ColorList[i] = lightgray;
    4220     else
    4221       ColorList[i] = white;
    4222   }
    4223   ShortestPathList[Root->nr] = 0;
    4224 
    4225   // and go on ... Queue always contains all lightgray vertices
    4226   while (!AtomStack->IsEmpty()) {
    4227     // we have to pop the oldest atom from stack. This keeps the atoms on the stack always of the same ShortestPath distance.
    4228     // e.g. if current atom is 2, push to end of stack are of length 3, but first all of length 2 would be popped. They again
    4229     // append length of 3 (their neighbours). Thus on stack we have always atoms of a certain length n at bottom of stack and
    4230     // followed by n+1 till top of stack.
    4231     Walker = AtomStack->PopFirst(); // pop oldest added
    4232     *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
    4233     for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
    4234       Binder = ListOfBondsPerAtom[Walker->nr][i];
    4235       if (Binder != NULL) { // don't look at bond equal NULL
    4236         OtherAtom = Binder->GetOtherAtom(Walker);
    4237         *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
    4238         if (ColorList[OtherAtom->nr] == white) {
    4239           if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
    4240             ColorList[OtherAtom->nr] = lightgray;
    4241           PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
    4242           ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
    4243           *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
    4244           if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
    4245             *out << Verbose(3);
    4246             if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
    4247               AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
    4248               *out << "Added OtherAtom " << OtherAtom->Name;
    4249               AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
    4250               AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
    4251               AddedBondList[Binder->nr]->Type = Binder->Type;
    4252               *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
    4253             } else {  // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
    4254               *out << "Not adding OtherAtom " << OtherAtom->Name;
    4255               if (AddedBondList[Binder->nr] == NULL) {
    4256                 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
    4257                 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
    4258                 AddedBondList[Binder->nr]->Type = Binder->Type;
    4259                 *out << ", added Bond " << *(AddedBondList[Binder->nr]);
    4260               } else
    4261                 *out << ", not added Bond ";
    4262             }
    4263             *out << ", putting OtherAtom into queue." << endl;
    4264             AtomStack->Push(OtherAtom);
    4265           } else { // out of bond order, then replace
    4266             if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
    4267               ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
    4268             if (Binder == Bond)
    4269               *out << Verbose(3) << "Not Queueing, is the Root bond";
    4270             else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
    4271               *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
    4272             if (!Binder->Cyclic)
    4273               *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
    4274             if (AddedBondList[Binder->nr] == NULL) {
    4275               if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
    4276                 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
    4277                 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
    4278                 AddedBondList[Binder->nr]->Type = Binder->Type;
    4279               } else {
    4280 #ifdef ADDHYDROGEN
    4281                 if (!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
    4282                   exit(1);
    4283 #endif
    4284               }
    4285             }
    4286           }
    4287         } else {
    4288           *out << Verbose(3) << "Not Adding, has already been visited." << endl;
    4289           // This has to be a cyclic bond, check whether it's present ...
    4290           if (AddedBondList[Binder->nr] == NULL) {
    4291             if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
    4292               AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
    4293               AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
    4294               AddedBondList[Binder->nr]->Type = Binder->Type;
    4295             } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
    4296 #ifdef ADDHYDROGEN
    4297               if(!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
    4298                 exit(1);
    4299 #endif
    4300             }
    4301           }
    4302         }
    4303       }
    4304     }
    4305     ColorList[Walker->nr] = black;
    4306     *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
    4307   }
    4308   Free(&PredecessorList);
    4309   Free(&ShortestPathList);
    4310   Free(&ColorList);
    4311   delete(AtomStack);
    4312 };
    4313 
    4314 /** Adds bond structure to this molecule from \a Father molecule.
    4315  * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
    4316  * with end points present in this molecule, bond is created in this molecule.
    4317  * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
    4318  * \param *out output stream for debugging
    4319  * \param *Father father molecule
    4320  * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
    4321  * \todo not checked, not fully working probably
    4322  */
    4323 bool molecule::BuildInducedSubgraph(ofstream *out, const molecule *Father)
    4324 {
    4325   atom *Walker = NULL, *OtherAtom = NULL;
    4326   bool status = true;
    4327   atom **ParentList = Malloc<atom*>(Father->AtomCount, "molecule::BuildInducedSubgraph: **ParentList");
    4328 
    4329   *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
    4330 
    4331   // reset parent list
    4332   *out << Verbose(3) << "Resetting ParentList." << endl;
    4333   for (int i=Father->AtomCount;i--;)
    4334     ParentList[i] = NULL;
    4335 
    4336   // fill parent list with sons
    4337   *out << Verbose(3) << "Filling Parent List." << endl;
    4338   Walker = start;
    4339   while (Walker->next != end) {
    4340     Walker = Walker->next;
    4341     ParentList[Walker->father->nr] = Walker;
    4342     // Outputting List for debugging
    4343     *out << Verbose(4) << "Son["<< Walker->father->nr <<"] of " << Walker->father <<  " is " << ParentList[Walker->father->nr] << "." << endl;
    4344   }
    4345 
    4346   // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
    4347   *out << Verbose(3) << "Creating bonds." << endl;
    4348   Walker = Father->start;
    4349   while (Walker->next != Father->end) {
    4350     Walker = Walker->next;
    4351     if (ParentList[Walker->nr] != NULL) {
    4352       if (ParentList[Walker->nr]->father != Walker) {
    4353         status = false;
    4354       } else {
    4355         for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
    4356           OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
    4357           if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
    4358             *out << Verbose(4) << "Endpoints of Bond " << Father->ListOfBondsPerAtom[Walker->nr][i] << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
    4359             AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], Father->ListOfBondsPerAtom[Walker->nr][i]->BondDegree);
    4360           }
    4361         }
    4362       }
    4363     }
    4364   }
    4365 
    4366   Free(&ParentList);
    4367   *out << Verbose(2) << "End of BuildInducedSubgraph." << endl;
    4368   return status;
    4369 };
    4370 
    4371 
    4372 /** Looks through a StackClass<atom *> and returns the likeliest removal candiate.
    4373  * \param *out output stream for debugging messages
    4374  * \param *&Leaf KeySet to look through
    4375  * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
    4376  * \param index of the atom suggested for removal
    4377  */
    4378 int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
    4379 {
    4380   atom *Runner = NULL;
    4381   int SP, Removal;
    4382 
    4383   *out << Verbose(2) << "Looking for removal candidate." << endl;
    4384   SP = -1; //0;  // not -1, so that Root is never removed
    4385   Removal = -1;
    4386   for (KeySet::iterator runner = Leaf->begin(); runner != Leaf->end(); runner++) {
    4387     Runner = FindAtom((*runner));
    4388     if (Runner->type->Z != 1) { // skip all those added hydrogens when re-filling snake stack
    4389       if (ShortestPathList[(*runner)] > SP) {  // remove the oldest one with longest shortest path
    4390         SP = ShortestPathList[(*runner)];
    4391         Removal = (*runner);
    4392       }
    4393     }
    4394   }
    4395   return Removal;
    4396 };
    4397 
    4398 /** Stores a fragment from \a KeySet into \a molecule.
    4399  * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
    4400  * molecule and adds missing hydrogen where bonds were cut.
    4401  * \param *out output stream for debugging messages
    4402  * \param &Leaflet pointer to KeySet structure
    4403  * \param IsAngstroem whether we have Ansgtroem or bohrradius
    4404  * \return pointer to constructed molecule
    4405  */
    4406 molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
    4407 {
    4408   atom *Runner = NULL, *FatherOfRunner = NULL, *OtherFather = NULL;
    4409   atom **SonList = Malloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList");
    4410   molecule *Leaf = new molecule(elemente);
    4411   bool LonelyFlag = false;
    4412   int size;
    4413 
    4414 //  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
    4415 
    4416   Leaf->BondDistance = BondDistance;
    4417   for(int i=NDIM*2;i--;)
    4418     Leaf->cell_size[i] = cell_size[i];
    4419 
    4420   // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
    4421   for(int i=AtomCount;i--;)
    4422     SonList[i] = NULL;
    4423 
    4424   // first create the minimal set of atoms from the KeySet
    4425   size = 0;
    4426   for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
    4427     FatherOfRunner = FindAtom((*runner));  // find the id
    4428     SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
    4429     size++;
    4430   }
    4431 
    4432   // create the bonds between all: Make it an induced subgraph and add hydrogen
    4433 //  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
    4434   Runner = Leaf->start;
    4435   while (Runner->next != Leaf->end) {
    4436     Runner = Runner->next;
    4437     LonelyFlag = true;
    4438     FatherOfRunner = Runner->father;
    4439     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
    4440       // create all bonds
    4441       for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
    4442         OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
    4443 //        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
    4444         if (SonList[OtherFather->nr] != NULL) {
    4445 //          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
    4446           if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
    4447 //            *out << Verbose(3) << "Adding Bond: ";
    4448 //            *out <<
    4449             Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
    4450 //            *out << "." << endl;
    4451             //NumBonds[Runner->nr]++;
    4452           } else {
    4453 //            *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
    4454           }
    4455           LonelyFlag = false;
    4456         } else {
    4457 //          *out << ", who has no son in this fragment molecule." << endl;
    4458 #ifdef ADDHYDROGEN
    4459           //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
    4460           if(!Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
    4461             exit(1);
    4462 #endif
    4463           //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
    4464         }
    4465       }
    4466     } else {
    4467       *out << Verbose(0) << "ERROR: Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
    4468     }
    4469     if ((LonelyFlag) && (size > 1)) {
    4470       *out << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl;
    4471     }
    4472 #ifdef ADDHYDROGEN
    4473     while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
    4474       Runner = Runner->next;
    4475 #endif
    4476   }
    4477   Leaf->CreateListOfBondsPerAtom(out);
    4478   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
    4479   Free(&SonList);
    4480 //  *out << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
    4481   return Leaf;
    4482 };
    4483 
    4484 /** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
    4485  * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
    4486  * computer game, that winds through the connected graph representing the molecule. Color (white,
    4487  * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
    4488  * creating only unique fragments and not additional ones with vertices simply in different sequence.
    4489  * The Predecessor is always the one that came before in discovering, needed on backstepping. And
    4490  * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
    4491  * stepping.
    4492  * \param *out output stream for debugging
    4493  * \param Order number of atoms in each fragment
    4494  * \param *configuration configuration for writing config files for each fragment
    4495  * \return List of all unique fragments with \a Order atoms
    4496  */
    4497 /*
    4498 MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
    4499 {
    4500   atom **PredecessorList = Malloc<atom*>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
    4501   int *ShortestPathList = Malloc<int>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
    4502   int *Labels = Malloc<int>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
    4503   enum Shading *ColorVertexList = Malloc<enum Shading>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
    4504   enum Shading *ColorEdgeList = Malloc<enum Shading>(BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
    4505   StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
    4506   StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
    4507   StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
    4508   MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
    4509   MoleculeListClass *FragmentList = NULL;
    4510   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
    4511   bond *Binder = NULL;
    4512   int RunningIndex = 0, FragmentCounter = 0;
    4513 
    4514   *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
    4515 
    4516   // reset parent list
    4517   *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
    4518   for (int i=0;i<AtomCount;i++) { // reset all atom labels
    4519     // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
    4520     Labels[i] = -1;
    4521     SonList[i] = NULL;
    4522     PredecessorList[i] = NULL;
    4523     ColorVertexList[i] = white;
    4524     ShortestPathList[i] = -1;
    4525   }
    4526   for (int i=0;i<BondCount;i++)
    4527     ColorEdgeList[i] = white;
    4528   RootStack->ClearStack();  // clearstack and push first atom if exists
    4529   TouchedStack->ClearStack();
    4530   Walker = start->next;
    4531   while ((Walker != end)
    4532 #ifdef ADDHYDROGEN
    4533    && (Walker->type->Z == 1)
    4534         }
    4535       }
    4536       *out << ", SP of " << ShortestPathList[Walker->nr]  << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
    4537 
    4538       // then check the stack for a newly stumbled upon fragment
    4539       if (SnakeStack->ItemCount() == Order) { // is stack full?
    4540         // store the fragment if it is one and get a removal candidate
    4541         Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
    4542         // remove the candidate if one was found
    4543         if (Removal != NULL) {
    4544           *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
    4545           SnakeStack->RemoveItem(Removal);
    4546           ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
    4547           if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
    4548             Walker = PredecessorList[Removal->nr];
    4549             *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
    4550           }
    4551         }
    4552       } else
    4553         Removal = NULL;
    4554 
    4555       // finally, look for a white neighbour as the next Walker
    4556       Binder = NULL;
    4557       if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {  // don't look, if a new walker has been set above
    4558         *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
    4559         OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
    4560         if (ShortestPathList[Walker->nr] < Order) {
    4561           for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
    4562             Binder = ListOfBondsPerAtom[Walker->nr][i];
    4563             *out << Verbose(2) << "Current bond is " << *Binder << ": ";
    4564             OtherAtom = Binder->GetOtherAtom(Walker);
    4565             if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
    4566               *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
    4567               //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
    4568             } else { // otherwise check its colour and element
    4569               if (
    4570               (OtherAtom->type->Z != 1) &&
    4571 #endif
    4572                     (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
    4573                 *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
    4574                 // i find it currently rather sensible to always set the predecessor in order to find one's way back
    4575                 //if (PredecessorList[OtherAtom->nr] == NULL) {
    4576                 PredecessorList[OtherAtom->nr] = Walker;
    4577                 *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
    4578                 //} else {
    4579                 //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
    4580                 //}
    4581                 Walker = OtherAtom;
    4582                 break;
    4583               } else {
    4584                 if (OtherAtom->type->Z == 1)
    4585                   *out << "Links to a hydrogen atom." << endl;
    4586                 else
    4587                   *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
    4588               }
    4589             }
    4590           }
    4591         } else {  // means we have stepped beyond the horizon: Return!
    4592           Walker = PredecessorList[Walker->nr];
    4593           OtherAtom = Walker;
    4594           *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
    4595         }
    4596         if (Walker != OtherAtom) {  // if no white neighbours anymore, color it black
    4597           *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
    4598           ColorVertexList[Walker->nr] = black;
    4599           Walker = PredecessorList[Walker->nr];
    4600         }
    4601       }
    4602     } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
    4603     *out << Verbose(2) << "Inner Looping is finished." << endl;
    4604 
    4605     // if we reset all AtomCount atoms, we have again technically O(N^2) ...
    4606     *out << Verbose(2) << "Resetting lists." << endl;
    4607     Walker = NULL;
    4608     Binder = NULL;
    4609     while (!TouchedStack->IsEmpty()) {
    4610       Walker = TouchedStack->PopLast();
    4611       *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
    4612       for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
    4613         ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
    4614       PredecessorList[Walker->nr] = NULL;
    4615       ColorVertexList[Walker->nr] = white;
    4616       ShortestPathList[Walker->nr] = -1;
    4617     }
    4618   }
    4619   *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
    4620 
    4621   // copy together
    4622   *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
    4623   FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
    4624   RunningIndex = 0;
    4625   while ((Leaflet != NULL) && (RunningIndex < FragmentCounter))  {
    4626     FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
    4627     Leaflet->Leaf = NULL; // prevent molecule from being removed
    4628     TempLeaf = Leaflet;
    4629     Leaflet = Leaflet->previous;
    4630     delete(TempLeaf);
    4631   };
    4632 
    4633   // free memory and exit
    4634   Free(&PredecessorList);
    4635   Free(&ShortestPathList);
    4636   Free(&Labels);
    4637   Free(&ColorVertexList);
    4638   delete(RootStack);
    4639   delete(TouchedStack);
    4640   delete(SnakeStack);
    4641 
    4642   *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
    4643   return FragmentList;
    4644 };
    4645 */
    4646 
    4647 /** Structure containing all values in power set combination generation.
    4648  */
    4649 struct UniqueFragments {
    4650   config *configuration;
    4651   atom *Root;
    4652   Graph *Leaflet;
    4653   KeySet *FragmentSet;
    4654   int ANOVAOrder;
    4655   int FragmentCounter;
    4656   int CurrentIndex;
    4657   double TEFactor;
    4658   int *ShortestPathList;
    4659   bool **UsedList;
    4660   bond **BondsPerSPList;
    4661   int *BondsPerSPCount;
    4662 };
    4663 
    4664 /** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.
    4665  * -# loops over every possible combination (2^dimension of edge set)
    4666  *  -# inserts current set, if there's still space left
    4667  *    -# yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root dist
    4668 ance+1
    4669  *    -# no: stores fragment into keyset list by calling InsertFragmentIntoGraph
    4670  *  -# removes all items added into the snake stack (in UniqueFragments structure) added during level (root
    4671 distance) and current set
    4672  * \param *out output stream for debugging
    4673  * \param FragmentSearch UniqueFragments structure with all values needed
    4674  * \param RootDistance current shortest path level, whose set of edges is represented by **BondsSet
    4675  * \param SetDimension Number of possible bonds on this level (i.e. size of the array BondsSet[])
    4676  * \param SubOrder remaining number of allowed vertices to add
    4677  */
    4678 void molecule::SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder)
    4679 {
    4680   atom *OtherWalker = NULL;
    4681   int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
    4682   int NumCombinations;
    4683   bool bit;
    4684   int bits, TouchedIndex, SubSetDimension, SP, Added;
    4685   int Removal;
    4686   int SpaceLeft;
    4687   int *TouchedList = Malloc<int>(SubOrder + 1, "molecule::SPFragmentGenerator: *TouchedList");
    4688   bond *Binder = NULL;
    4689   bond **BondsList = NULL;
    4690   KeySetTestPair TestKeySetInsert;
    4691 
    4692   NumCombinations = 1 << SetDimension;
    4693 
    4694   // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
    4695   // von Endstuecken (aus den Bonds) hinzugefuegt werden und fuer verbleibende ANOVAOrder
    4696   // rekursiv GraphCrawler in der naechsten Ebene aufgerufen werden
    4697 
    4698   *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
    4699   *out << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->Root << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " <<  NumCombinations-1 << " combination(s)." << endl;
    4700 
    4701   // initialised touched list (stores added atoms on this level)
    4702   *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
    4703   for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
    4704     TouchedList[TouchedIndex] = -1;
    4705   TouchedIndex = 0;
    4706 
    4707   // create every possible combination of the endpieces
    4708   *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
    4709   for (int i=1;i<NumCombinations;i++) {  // sweep through all power set combinations (skip empty set!)
    4710     // count the set bit of i
    4711     bits = 0;
    4712     for (int j=SetDimension;j--;)
    4713       bits += (i & (1 << j)) >> j;
    4714 
    4715     *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
    4716     if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
    4717       // --1-- add this set of the power set of bond partners to the snake stack
    4718       Added = 0;
    4719       for (int j=0;j<SetDimension;j++) {  // pull out every bit by shifting
    4720         bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
    4721         if (bit) {  // if bit is set, we add this bond partner
    4722           OtherWalker = BondsSet[j]->rightatom;  // rightatom is always the one more distant, i.e. the one to add
    4723           //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
    4724           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
    4725           TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
    4726           if (TestKeySetInsert.second) {
    4727             TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
    4728             Added++;
    4729           } else {
    4730             *out << Verbose(2+verbosity) << "This was item was already present in the keyset." << endl;
    4731           }
    4732             //FragmentSearch->UsedList[OtherWalker->nr][i] = true;
    4733           //}
    4734         } else {
    4735           *out << Verbose(2+verbosity) << "Not adding." << endl;
    4736         }
    4737       }
    4738 
    4739       SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
    4740       if (SpaceLeft > 0) {
    4741         *out << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
    4742         if (SubOrder > 1) {    // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
    4743           // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
    4744           SP = RootDistance+1;  // this is the next level
    4745           // first count the members in the subset
    4746           SubSetDimension = 0;
    4747           Binder = FragmentSearch->BondsPerSPList[2*SP];    // start node for this level
    4748           while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {    // compare to end node of this level
    4749             Binder = Binder->next;
    4750             for (int k=TouchedIndex;k--;) {
    4751               if (Binder->Contains(TouchedList[k]))   // if we added this very endpiece
    4752                 SubSetDimension++;
    4753             }
    4754           }
    4755           // then allocate and fill the list
    4756           BondsList = Malloc<bond*>(SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
    4757           SubSetDimension = 0;
    4758           Binder = FragmentSearch->BondsPerSPList[2*SP];
    4759           while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {
    4760             Binder = Binder->next;
    4761             for (int k=0;k<TouchedIndex;k++) {
    4762               if (Binder->leftatom->nr == TouchedList[k])   // leftatom is always the close one
    4763                 BondsList[SubSetDimension++] = Binder;
    4764             }
    4765           }
    4766           *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
    4767           SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
    4768           Free(&BondsList);
    4769         }
    4770       } else {
    4771         // --2-- otherwise store the complete fragment
    4772         *out << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
    4773         // store fragment as a KeySet
    4774         *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
    4775         for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
    4776           *out << (*runner) << " ";
    4777         *out << endl;
    4778         //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
    4779           //*out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl;
    4780         InsertFragmentIntoGraph(out, FragmentSearch);
    4781         //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
    4782         //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
    4783       }
    4784 
    4785       // --3-- remove all added items in this level from snake stack
    4786       *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
    4787       for(int j=0;j<TouchedIndex;j++) {
    4788         Removal = TouchedList[j];
    4789         *out << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
    4790         FragmentSearch->FragmentSet->erase(Removal);
    4791         TouchedList[j] = -1;
    4792       }
    4793       *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
    4794       for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
    4795         *out << (*runner) << " ";
    4796       *out << endl;
    4797       TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
    4798     } else {
    4799       *out << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
    4800     }
    4801   }
    4802   Free(&TouchedList);
    4803   *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
    4804 };
    4805 
    4806 /** For a given keyset \a *Fragment, checks whether it is connected in the current molecule.
    4807  * \param *out output stream for debugging
    4808  * \param *Fragment Keyset of fragment's vertices
    4809  * \return true - connected, false - disconnected
    4810  * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
    4811  */
    4812 bool molecule::CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment)
    4813 {
    4814   atom *Walker = NULL, *Walker2 = NULL;
    4815   bool BondStatus = false;
    4816   int size;
    4817 
    4818   *out << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl;
    4819   *out << Verbose(2) << "Disconnected atom: ";
    4820 
    4821   // count number of atoms in graph
    4822   size = 0;
    4823   for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
    4824     size++;
    4825   if (size > 1)
    4826     for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++) {
    4827       Walker = FindAtom(*runner);
    4828       BondStatus = false;
    4829       for(KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
    4830         Walker2 = FindAtom(*runners);
    4831         for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
    4832           if (ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker) == Walker2) {
    4833             BondStatus = true;
    4834             break;
    4835           }
    4836           if (BondStatus)
    4837             break;
    4838         }
    4839       }
    4840       if (!BondStatus) {
    4841         *out << (*Walker) << endl;
    4842         return false;
    4843       }
    4844     }
    4845   else {
    4846     *out << "none." << endl;
    4847     return true;
    4848   }
    4849   *out << "none." << endl;
    4850 
    4851   *out << Verbose(1) << "End of CheckForConnectedSubgraph" << endl;
    4852 
    4853   return true;
    4854 }
    4855 
    4856 /** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
    4857  * -# initialises UniqueFragments structure
    4858  * -# fills edge list via BFS
    4859  * -# creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as
    4860  root distance, the edge set, its dimension and the current suborder
    4861  * -# Free'ing structure
    4862  * Note that we may use the fact that the atoms are SP-ordered on the atomstack. I.e. when popping always the last, we first get all
    4863  * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
    4864  * \param *out output stream for debugging
    4865  * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
    4866  * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
    4867  * \param RestrictedKeySet Restricted vertex set to use in context of molecule
    4868  * \return number of inserted fragments
    4869  * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
    4870  */
    4871 int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
    4872 {
    4873   int SP, AtomKeyNr;
    4874   atom *Walker = NULL, *OtherWalker = NULL, *Predecessor = NULL;
    4875   bond *Binder = NULL;
    4876   bond *CurrentEdge = NULL;
    4877   bond **BondsList = NULL;
    4878   int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
    4879   int Counter = FragmentSearch.FragmentCounter;
    4880   int RemainingWalkers;
    4881 
    4882   *out << endl;
    4883   *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;
    4884 
    4885   // prepare Label and SP arrays of the BFS search
    4886   FragmentSearch.ShortestPathList[FragmentSearch.Root->nr] = 0;
    4887 
    4888   // prepare root level (SP = 0) and a loop bond denoting Root
    4889   for (int i=1;i<Order;i++)
    4890     FragmentSearch.BondsPerSPCount[i] = 0;
    4891   FragmentSearch.BondsPerSPCount[0] = 1;
    4892   Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
    4893   add(Binder, FragmentSearch.BondsPerSPList[1]);
    4894 
    4895   // do a BFS search to fill the SP lists and label the found vertices
    4896   // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
    4897   // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
    4898   // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
    4899   // (EdgeinSPLevel) of this tree ...
    4900   // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
    4901   // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
    4902   *out << endl;
    4903   *out << Verbose(0) << "Starting BFS analysis ..." << endl;
    4904   for (SP = 0; SP < (Order-1); SP++) {
    4905     *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with " << FragmentSearch.BondsPerSPCount[SP] << " item(s)";
    4906     if (SP > 0) {
    4907       *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
    4908       FragmentSearch.BondsPerSPCount[SP] = 0;
    4909     } else
    4910       *out << "." << endl;
    4911 
    4912     RemainingWalkers = FragmentSearch.BondsPerSPCount[SP];
    4913     CurrentEdge = FragmentSearch.BondsPerSPList[2*SP];    /// start of this SP level's list
    4914     while (CurrentEdge->next != FragmentSearch.BondsPerSPList[2*SP+1]) {    /// end of this SP level's list
    4915       CurrentEdge = CurrentEdge->next;
    4916       RemainingWalkers--;
    4917       Walker = CurrentEdge->rightatom;    // rightatom is always the one more distant
    4918       Predecessor = CurrentEdge->leftatom;    // ... and leftatom is predecessor
    4919       AtomKeyNr = Walker->nr;
    4920       *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
    4921       // check for new sp level
    4922       // go through all its bonds
    4923       *out << Verbose(1) << "Going through all bonds of Walker." << endl;
    4924       for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
    4925         Binder = ListOfBondsPerAtom[AtomKeyNr][i];
    4926         OtherWalker = Binder->GetOtherAtom(Walker);
    4927         if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
    4928   #ifdef ADDHYDROGEN
    4929          && (OtherWalker->type->Z != 1)
    4930   #endif
    4931                                                               ) {  // skip hydrogens and restrict to fragment
    4932           *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
    4933           // set the label if not set (and push on root stack as well)
    4934           if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
    4935             FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
    4936             *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
    4937             // add the bond in between to the SP list
    4938             Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
    4939             add(Binder, FragmentSearch.BondsPerSPList[2*(SP+1)+1]);
    4940             FragmentSearch.BondsPerSPCount[SP+1]++;
    4941             *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl;
    4942           } else {
    4943             if (OtherWalker != Predecessor)
    4944               *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->nr << " is smaller than that of Root " << RootKeyNr << "." << endl;
    4945             else
    4946               *out << Verbose(3) << "This is my predecessor " << *Predecessor << "." << endl;
    4947           }
    4948         } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
    4949       }
    4950     }
    4951   }
    4952 
    4953   // outputting all list for debugging
    4954   *out << Verbose(0) << "Printing all found lists." << endl;
    4955   for(int i=1;i<Order;i++) {    // skip the root edge in the printing
    4956     Binder = FragmentSearch.BondsPerSPList[2*i];
    4957     *out << Verbose(1) << "Current SP level is " << i << "." << endl;
    4958     while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
    4959       Binder = Binder->next;
    4960       *out << Verbose(2) << *Binder << endl;
    4961     }
    4962   }
    4963 
    4964   // creating fragments with the found edge sets  (may be done in reverse order, faster)
    4965   SP = -1;  // the Root <-> Root edge must be subtracted!
    4966   for(int i=Order;i--;) { // sum up all found edges
    4967     Binder = FragmentSearch.BondsPerSPList[2*i];
    4968     while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
    4969       Binder = Binder->next;
    4970       SP ++;
    4971     }
    4972   }
    4973   *out << Verbose(0) << "Total number of edges is " << SP << "." << endl;
    4974   if (SP >= (Order-1)) {
    4975     // start with root (push on fragment stack)
    4976     *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
    4977     FragmentSearch.FragmentSet->clear();
    4978     *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
    4979     // prepare the subset and call the generator
    4980     BondsList = Malloc<bond*>(FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
    4981     BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond
    4982 
    4983     SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
    4984 
    4985     Free(&BondsList);
    4986   } else {
    4987     *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
    4988   }
    4989 
    4990   // as FragmentSearch structure is used only once, we don't have to clean it anymore
    4991   // remove root from stack
    4992   *out << Verbose(0) << "Removing root again from stack." << endl;
    4993   FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
    4994 
    4995   // free'ing the bonds lists
    4996   *out << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl;
    4997   for(int i=Order;i--;) {
    4998     *out << Verbose(1) << "Current SP level is " << i << ": ";
    4999     Binder = FragmentSearch.BondsPerSPList[2*i];
    5000     while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
    5001       Binder = Binder->next;
    5002       // *out << "Removing atom " << Binder->leftatom->nr << " and " << Binder->rightatom->nr << "." << endl; // make sure numbers are local
    5003       FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
    5004       FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
    5005     }
    5006     // delete added bonds
    5007     cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
    5008     // also start and end node
    5009     *out << "cleaned." << endl;
    5010   }
    5011 
    5012   // return list
    5013   *out << Verbose(0) << "End of PowerSetGenerator." << endl;
    5014   return (FragmentSearch.FragmentCounter - Counter);
    5015 };
    5016 
    5017 /** Corrects the nuclei position if the fragment was created over the cell borders.
    5018  * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
    5019  * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
    5020  * and re-add the bond. Looping on the distance check.
    5021  * \param *out ofstream for debugging messages
    5022  */
    5023 void molecule::ScanForPeriodicCorrection(ofstream *out)
    5024 {
    5025   bond *Binder = NULL;
    5026   bond *OtherBinder = NULL;
    5027   atom *Walker = NULL;
    5028   atom *OtherWalker = NULL;
    5029   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    5030   enum Shading *ColorList = NULL;
    5031   double tmp;
    5032   Vector Translationvector;
    5033   //class StackClass<atom *> *CompStack = NULL;
    5034   class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
    5035   bool flag = true;
    5036 
    5037   *out << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl;
    5038 
    5039   ColorList = Malloc<enum Shading>(AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");
    5040   while (flag) {
    5041     // remove bonds that are beyond bonddistance
    5042     for(int i=NDIM;i--;)
    5043       Translationvector.x[i] = 0.;
    5044     // scan all bonds
    5045     Binder = first;
    5046     flag = false;
    5047     while ((!flag) && (Binder->next != last)) {
    5048       Binder = Binder->next;
    5049       for (int i=NDIM;i--;) {
    5050         tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
    5051         //*out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
    5052         if (tmp > BondDistance) {
    5053           OtherBinder = Binder->next; // note down binding partner for later re-insertion
    5054           unlink(Binder);   // unlink bond
    5055           *out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
    5056           flag = true;
    5057           break;
    5058         }
    5059       }
    5060     }
    5061     if (flag) {
    5062       // create translation vector from their periodically modified distance
    5063       for (int i=NDIM;i--;) {
    5064         tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
    5065         if (fabs(tmp) > BondDistance)
    5066           Translationvector.x[i] = (tmp < 0) ? +1. : -1.;
    5067       }
    5068       Translationvector.MatrixMultiplication(matrix);
    5069       //*out << Verbose(3) << "Translation vector is ";
    5070       Translationvector.Output(out);
    5071       *out << endl;
    5072       // apply to all atoms of first component via BFS
    5073       for (int i=AtomCount;i--;)
    5074         ColorList[i] = white;
    5075       AtomStack->Push(Binder->leftatom);
    5076       while (!AtomStack->IsEmpty()) {
    5077         Walker = AtomStack->PopFirst();
    5078         //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
    5079         ColorList[Walker->nr] = black;    // mark as explored
    5080         Walker->x.AddVector(&Translationvector); // translate
    5081         for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {  // go through all binding partners
    5082           if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
    5083             OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
    5084             if (ColorList[OtherWalker->nr] == white) {
    5085               AtomStack->Push(OtherWalker); // push if yet unexplored
    5086             }
    5087           }
    5088         }
    5089       }
    5090       // re-add bond
    5091       link(Binder, OtherBinder);
    5092     } else {
    5093       *out << Verbose(3) << "No corrections for this fragment." << endl;
    5094     }
    5095     //delete(CompStack);
    5096   }
    5097 
    5098   // free allocated space from ReturnFullMatrixforSymmetric()
    5099   delete(AtomStack);
    5100   Free(&ColorList);
    5101   Free(&matrix);
    5102   *out << Verbose(2) << "End of ScanForPeriodicCorrection." << endl;
    5103 };
    51041041
    51051042/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
    … …  
    51221059};
    51231060
    5124 bool KeyCompare::operator() (const KeySet SubgraphA, const KeySet SubgraphB) const
    5125 {
    5126   //cout << "my check is used." << endl;
    5127   if (SubgraphA.size() < SubgraphB.size()) {
    5128     return true;
    5129   } else {
    5130     if (SubgraphA.size() > SubgraphB.size()) {
    5131       return false;
    5132     } else {
    5133       KeySet::iterator IteratorA = SubgraphA.begin();
    5134       KeySet::iterator IteratorB = SubgraphB.begin();
    5135       while ((IteratorA != SubgraphA.end()) && (IteratorB != SubgraphB.end())) {
    5136         if ((*IteratorA) <  (*IteratorB))
    5137           return true;
    5138         else if ((*IteratorA) > (*IteratorB)) {
    5139             return false;
    5140           } // else, go on to next index
    5141         IteratorA++;
    5142         IteratorB++;
    5143       } // end of while loop
    5144     }// end of check in case of equal sizes
    5145   }
    5146   return false; // if we reach this point, they are equal
    5147 };
    5148 
    5149 //bool operator < (KeySet SubgraphA, KeySet SubgraphB)
    5150 //{
    5151 //  return KeyCompare(SubgraphA, SubgraphB);
    5152 //};
    5153 
    5154 /** Checking whether KeySet is not already present in Graph, if so just adds factor.
    5155  * \param *out output stream for debugging
    5156  * \param &set KeySet to insert
    5157  * \param &graph Graph to insert into
    5158  * \param *counter pointer to unique fragment count
    5159  * \param factor energy factor for the fragment
    5160  */
    5161 inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment)
    5162 {
    5163   GraphTestPair testGraphInsert;
    5164 
    5165   testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,Fragment->TEFactor)));  // store fragment number and current factor
    5166   if (testGraphInsert.second) {
    5167     *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " successfully inserted." << endl;
    5168     Fragment->FragmentCounter++;
    5169   } else {
    5170     *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
    5171     ((*(testGraphInsert.first)).second).second += Fragment->TEFactor;  // increase the "created" counter
    5172     *out << Verbose(2) << "New factor is " << ((*(testGraphInsert.first)).second).second << "." << endl;
    5173   }
    5174 };
    5175 //void inline InsertIntoGraph(ofstream *out, KeyStack &stack, Graph &graph, int *counter, double factor)
    5176 //{
    5177 //  // copy stack contents to set and call overloaded function again
    5178 //  KeySet set;
    5179 //  for(KeyStack::iterator runner = stack.begin(); runner != stack.begin(); runner++)
    5180 //    set.insert((*runner));
    5181 //  InsertIntoGraph(out, set, graph, counter, factor);
    5182 //};
    5183 
    5184 /** Inserts each KeySet in \a graph2 into \a graph1.
    5185  * \param *out output stream for debugging
    5186  * \param graph1 first (dest) graph
    5187  * \param graph2 second (source) graph
    5188  * \param *counter keyset counter that gets increased
    5189  */
    5190 inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
    5191 {
    5192   GraphTestPair testGraphInsert;
    5193 
    5194   for(Graph::iterator runner = graph2.begin(); runner != graph2.end(); runner++) {
    5195     testGraphInsert = graph1.insert(GraphPair ((*runner).first,pair<int,double>((*counter)++,((*runner).second).second)));  // store fragment number and current factor
    5196     if (testGraphInsert.second) {
    5197       *out << Verbose(2) << "KeySet " << (*counter)-1 << " successfully inserted." << endl;
    5198     } else {
    5199       *out << Verbose(2) << "KeySet " << (*counter)-1 << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
    5200       ((*(testGraphInsert.first)).second).second += (*runner).second.second;
    5201       *out << Verbose(2) << "New factor is " << (*(testGraphInsert.first)).second.second << "." << endl;
    5202     }
    5203   }
    5204 };
    5205 
    5206 
    5207 /** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
    5208  * -# constructs a complete keyset of the molecule
    5209  * -# In a loop over all possible roots from the given rootstack
    5210  *  -# increases order of root site
    5211  *  -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
    5212  *  -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
    5213 as the restricted one and each site in the set as the root)
    5214  *  -# these are merged into a fragment list of keysets
    5215  * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
    5216  * Important only is that we create all fragments, it is not important if we create them more than once
    5217  * as these copies are filtered out via use of the hash table (KeySet).
    5218  * \param *out output stream for debugging
    5219  * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
    5220  * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
    5221  * \param *MinimumRingSize minimum ring size for each atom (molecule::Atomcount)
    5222  * \return pointer to Graph list
    5223  */
    5224 void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
    5225 {
    5226   Graph ***FragmentLowerOrdersList = NULL;
    5227   int NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
    5228   int counter = 0, Order;
    5229   int UpgradeCount = RootStack.size();
    5230   KeyStack FragmentRootStack;
    5231   int RootKeyNr, RootNr;
    5232   struct UniqueFragments FragmentSearch;
    5233 
    5234   *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
    5235 
    5236   // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
    5237   // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
    5238   NumMoleculesOfOrder = Malloc<int>(UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
    5239   FragmentLowerOrdersList = Malloc<Graph**>(UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
    5240 
    5241   // initialise the fragments structure
    5242   FragmentSearch.ShortestPathList = Malloc<int>(AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
    5243   FragmentSearch.FragmentCounter = 0;
    5244   FragmentSearch.FragmentSet = new KeySet;
    5245   FragmentSearch.Root = FindAtom(RootKeyNr);
    5246   for (int i=AtomCount;i--;) {
    5247     FragmentSearch.ShortestPathList[i] = -1;
    5248   }
    5249 
    5250   // Construct the complete KeySet which we need for topmost level only (but for all Roots)
    5251   atom *Walker = start;
    5252   KeySet CompleteMolecule;
    5253   while (Walker->next != end) {
    5254     Walker = Walker->next;
    5255     CompleteMolecule.insert(Walker->GetTrueFather()->nr);
    5256   }
    5257 
    5258   // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
    5259   // each has to be split up. E.g. for the second level we have one from 5th, one from 4th, two from 3th (which in turn is one from 5th, one from 4th),
    5260   // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
    5261   // with the order along the cells as this: 5433222211111111 for BondOrder 5 needing 16=pow(2,5-1) cells (only we use bit-shifting which is faster)
    5262   RootNr = 0;   // counts through the roots in RootStack
    5263   while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
    5264     RootKeyNr = RootStack.front();
    5265     RootStack.pop_front();
    5266     Walker = FindAtom(RootKeyNr);
    5267     // check cyclic lengths
    5268     //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
    5269     //  *out << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
    5270     //} else
    5271     {
    5272       // increase adaptive order by one
    5273       Walker->GetTrueFather()->AdaptiveOrder++;
    5274       Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
    5275 
    5276       // initialise Order-dependent entries of UniqueFragments structure
    5277       FragmentSearch.BondsPerSPList = Malloc<bond*>(Order * 2, "molecule::PowerSetGenerator: ***BondsPerSPList");
    5278       FragmentSearch.BondsPerSPCount = Malloc<int>(Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
    5279       for (int i=Order;i--;) {
    5280         FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
    5281         FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
    5282         FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
    5283         FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
    5284         FragmentSearch.BondsPerSPCount[i] = 0;
    5285       }
    5286 
    5287       // allocate memory for all lower level orders in this 1D-array of ptrs
    5288       NumLevels = 1 << (Order-1); // (int)pow(2,Order);
    5289       FragmentLowerOrdersList[RootNr] = Malloc<Graph*>(NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
    5290       for (int i=0;i<NumLevels;i++)
    5291         FragmentLowerOrdersList[RootNr][i] = NULL;
    5292 
    5293       // create top order where nothing is reduced
    5294       *out << Verbose(0) << "==============================================================================================================" << endl;
    5295       *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl; // , NumLevels is " << NumLevels << "
    5296 
    5297       // Create list of Graphs of current Bond Order (i.e. F_{ij})
    5298       FragmentLowerOrdersList[RootNr][0] =  new Graph;
    5299       FragmentSearch.TEFactor = 1.;
    5300       FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0];      // set to insertion graph
    5301       FragmentSearch.Root = Walker;
    5302       NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
    5303       *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    5304       if (NumMoleculesOfOrder[RootNr] != 0) {
    5305         NumMolecules = 0;
    5306 
    5307         // we don't have to dive into suborders! These keysets are all already created on lower orders!
    5308         // this was all ancient stuff, when we still depended on the TEFactors (and for those the suborders were needed)
    5309 
    5310 //        if ((NumLevels >> 1) > 0) {
    5311 //          // create lower order fragments
    5312 //          *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
    5313 //          Order = Walker->AdaptiveOrder;
    5314 //          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
    5315 //            // step down to next order at (virtual) boundary of powers of 2 in array
    5316 //            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
    5317 //              Order--;
    5318 //            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
    5319 //            for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
    5320 //              int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
    5321 //              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
    5322 //              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
    5323 //
    5324 //              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
    5325 //              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
    5326 //              //NumMolecules = 0;
    5327 //              FragmentLowerOrdersList[RootNr][dest] = new Graph;
    5328 //              for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
    5329 //                for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    5330 //                  Graph TempFragmentList;
    5331 //                  FragmentSearch.TEFactor = -(*runner).second.second;
    5332 //                  FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
    5333 //                  FragmentSearch.Root = FindAtom(*sprinter);
    5334 //                  NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
    5335 //                  // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
    5336 //                  *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
    5337 //                  InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
    5338 //                }
    5339 //              }
    5340 //              *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
    5341 //            }
    5342 //          }
    5343 //        }
    5344       } else {
    5345         Walker->GetTrueFather()->MaxOrder = true;
    5346 //        *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
    5347       }
    5348       // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
    5349       //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
    5350       TotalNumMolecules += NumMoleculesOfOrder[RootNr];
    5351 //      *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    5352       RootStack.push_back(RootKeyNr); // put back on stack
    5353       RootNr++;
    5354 
    5355       // free Order-dependent entries of UniqueFragments structure for next loop cycle
    5356       Free(&FragmentSearch.BondsPerSPCount);
    5357       for (int i=Order;i--;) {
    5358         delete(FragmentSearch.BondsPerSPList[2*i]);
    5359         delete(FragmentSearch.BondsPerSPList[2*i+1]);
    5360       }
    5361       Free(&FragmentSearch.BondsPerSPList);
    5362     }
    5363   }
    5364   *out << Verbose(0) << "==============================================================================================================" << endl;
    5365   *out << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl;
    5366   *out << Verbose(0) << "==============================================================================================================" << endl;
    5367 
    5368   // cleanup FragmentSearch structure
    5369   Free(&FragmentSearch.ShortestPathList);
    5370   delete(FragmentSearch.FragmentSet);
    5371 
    5372   // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
    5373   // 5433222211111111
    5374   // 43221111
    5375   // 3211
    5376   // 21
    5377   // 1
    5378 
    5379   // Subsequently, we combine all into a single list (FragmentList)
    5380 
    5381   *out << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl;
    5382   if (FragmentList == NULL) {
    5383     FragmentList = new Graph;
    5384     counter = 0;
    5385   } else {
    5386     counter = FragmentList->size();
    5387   }
    5388   RootNr = 0;
    5389   while (!RootStack.empty()) {
    5390     RootKeyNr = RootStack.front();
    5391     RootStack.pop_front();
    5392     Walker = FindAtom(RootKeyNr);
    5393     NumLevels = 1 << (Walker->AdaptiveOrder - 1);
    5394     for(int i=0;i<NumLevels;i++) {
    5395       if (FragmentLowerOrdersList[RootNr][i] != NULL) {
    5396         InsertGraphIntoGraph(out, *FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
    5397         delete(FragmentLowerOrdersList[RootNr][i]);
    5398       }
    5399     }
    5400     Free(&FragmentLowerOrdersList[RootNr]);
    5401     RootNr++;
    5402   }
    5403   Free(&FragmentLowerOrdersList);
    5404   Free(&NumMoleculesOfOrder);
    5405 
    5406   *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
    5407 };
    54081061
    54091062/** Comparison function for GSL heapsort on distances in two molecules.

  • src/molecule.hpp

    rd09ff7 rcee0b57  
    1 /** \file molecules.hpp
     1/** \file molecule.hpp
    22 *
    33 * Class definitions of atom and molecule, element and periodentafel
    … …  
    133133
    134134  // templates for allowing global manipulation of all vectors
     135  template <typename res> void ActOnAllVectors( res (Vector::*f)() );
    135136  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t );
    136137  template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u );
    137138  template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v);
     139
     140  // templates for allowing global manipulation of all atoms
     141  template <typename res> void ActOnAllAtoms( res (molecule::*f)(atom *) );
     142  template <typename res> void ActOnAllAtoms( res (atom::*f)() );
     143  template <typename res, typename T> void ActOnAllAtoms( res (atom::*f)(T), T t );
     144  template <typename res, typename T, typename U> void ActOnAllAtoms( res (atom::*f)(T, U), T t, U u );
     145  template <typename res, typename T, typename U, typename V> void ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v);
    138146
    139147  /// remove atoms from molecule.
    … …  
    147155  bool AddXYZFile(string filename);
    148156  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
    149   bond * AddBond(atom *first, atom *second, int degree);
     157  bond * AddBond(atom *first, atom *second, int degree = 1);
    150158  bool RemoveBond(bond *pointer);
    151159  bool RemoveBonds(atom *BondPartner);
    … …  
    255263
    256264
     265template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() ) {
     266  atom *Walker = start;
     267  while (Walker->next != end) {
     268    Walker = Walker->next;
     269    ((Walker->node)->*f)();
     270  }
     271};
    257272template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t )
    258273{
    … …  
    260275  while (Walker->next != end) {
    261276    Walker = Walker->next;
    262     ((*Walker->node)->*f)(t);
    263   }
    264 };
    265 
     277    ((Walker->node)->*f)(t);
     278  }
     279};
    266280template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u )
    267281{
    … …  
    269283  while (Walker->next != end) {
    270284    Walker = Walker->next;
    271     ((*Walker->node)->*f)(t, u);
    272   }
    273 };
    274 
     285    ((Walker->node)->*f)(t, u);
     286  }
     287};
    275288template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v)
    276289{
    … …  
    278291  while (Walker->next != end) {
    279292    Walker = Walker->next;
    280     ((*Walker->node)->*f)(t, u, v);
     293    ((Walker->node)->*f)(t, u, v);
     294  }
     295};
     296
     297template <typename res> void molecule::ActOnAllAtoms( res (molecule::*f)(atom *)) {
     298  atom *Walker = start;
     299  while (Walker->next != end) {
     300    Walker = Walker->next;
     301    (*f)(Walker);
     302  }
     303};
     304
     305template <typename res> void molecule::ActOnAllAtoms( res (atom::*f)()) {
     306  atom *Walker = start;
     307  while (Walker->next != end) {
     308    Walker = Walker->next;
     309    (Walker->*f)();
     310  }
     311};
     312template <typename res, typename T> void molecule::ActOnAllAtoms( res (atom::*f)(T), T t )
     313{
     314  atom *Walker = start;
     315  while (Walker->next != end) {
     316    Walker = Walker->next;
     317    (Walker->*f)(t);
     318  }
     319};
     320template <typename res, typename T, typename U> void molecule::ActOnAllAtoms( res (atom::*f)(T, U), T t, U u )
     321{
     322  atom *Walker = start;
     323  while (Walker->next != end) {
     324    Walker = Walker->next;
     325    (Walker->*f)(t, u);
     326  }
     327};
     328template <typename res, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (atom::*f)(T, U, V), T t, U u, V v)
     329{
     330  atom *Walker = start;
     331  while (Walker->next != end) {
     332    Walker = Walker->next;
     333    (Walker->*f)(t, u, v);
    281334  }
    282335};

  • src/moleculelist.cpp

    rd09ff7 rcee0b57  
    77#include "boundary.hpp"
    88#include "config.hpp"
    9 #include "molecules.hpp"
     9#include "molecule.hpp"
    1010#include "memoryallocator.hpp"
    1111

  • src/verbose.cpp

    rd09ff7 rcee0b57  
    1 #include "molecules.hpp"
     1#include "molecule.hpp"
    22
    33/** Prints the tabs according to verbosity stored in the temporary constructed class.
Note: See TracChangeset for help on using the changeset viewer.