Changeset 69eb71 for src/molecules.cpp


Ignore:
Timestamp:
Dec 3, 2008, 2:12:05 PM (16 years ago)
Author:
Christian Neuen <neuen@…>
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:
f683fe
Parents:
1ffa21
Message:

Multiple changes to boundary, currently not fully operational.
Molecules has a new routine to create adjacency lists, reading bonds from a dbond file instead of looking for the distances by itself.
Vector function Project onto plane has been updated.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecules.cpp

    r1ffa21 r69eb71  
    11/** \file molecules.cpp
    2  * 
     2 *
    33 * Functions for the class molecule.
    4  * 
     4 *
    55 */
    66
    … …  
    2525      sum += (gsl_vector_get(x,j) - (vectors[i])->x[j])*(gsl_vector_get(x,j) - (vectors[i])->x[j]);
    2626  }
    27  
     27
    2828  return sum;
    2929};
    … …  
    3434 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
    3535 */
    36 molecule::molecule(periodentafel *teil) 
    37 { 
     36molecule::molecule(periodentafel *teil)
     37{
    3838  // init atom chain list
    39   start = new atom; 
     39  start = new atom;
    4040  end = new atom;
    41   start->father = NULL; 
     41  start->father = NULL;
    4242  end->father = NULL;
    4343  link(start,end);
    … …  
    4646  last = new bond(start, end, 1, -1);
    4747  link(first,last);
    48   // other stuff 
     48  // other stuff
    4949  MDSteps = 0;
    50   last_atom = 0; 
     50  last_atom = 0;
    5151  elemente = teil;
    5252  AtomCount = 0;
    … …  
    6767 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
    6868 */
    69 molecule::~molecule() 
     69molecule::~molecule()
    7070{
    7171  if (ListOfBondsPerAtom != NULL)
    … …  
    7878  delete(last);
    7979  delete(end);
    80   delete(start); 
     80  delete(start);
    8181};
    8282
    8383/** Adds given atom \a *pointer from molecule list.
    84  * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount 
     84 * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
    8585 * \param *pointer allocated and set atom
    8686 * \return true - succeeded, false - atom not found in list
    8787 */
    8888bool molecule::AddAtom(atom *pointer)
    89 { 
     89{
    9090  if (pointer != NULL) {
    91     pointer->sort = &pointer->nr; 
     91    pointer->sort = &pointer->nr;
    9292    pointer->nr = last_atom++;  // increase number within molecule
    9393    AtomCount++;
    … …  
    106106    return add(pointer, end);
    107107  } else
    108     return false; 
     108    return false;
    109109};
    110110
    … …  
    115115 */
    116116atom * molecule::AddCopyAtom(atom *pointer)
    117 { 
     117{
    118118  if (pointer != NULL) {
    119119        atom *walker = new atom();
    … …  
    122122    walker->v.CopyVector(&pointer->v); // copy velocity
    123123    walker->FixedIon = pointer->FixedIon;
    124     walker->sort = &walker->nr; 
     124    walker->sort = &walker->nr;
    125125    walker->nr = last_atom++;  // increase number within molecule
    126126    walker->father = pointer; //->GetTrueFather();
    … …  
    133133    return walker;
    134134  } else
    135     return NULL; 
     135    return NULL;
    136136};
    137137
    … …  
    156156 *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
    157157 *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
    158  *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2} 
     158 *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
    159159 *    \f]
    160  *    as the coordination of all three atoms in the coordinate system of these three vectors: 
     160 *    as the coordination of all three atoms in the coordinate system of these three vectors:
    161161 *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
    162  * 
     162 *
    163163 * \param *out output stream for debugging
    164  * \param *Bond pointer to bond between \a *origin and \a *replacement 
    165  * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin) 
     164 * \param *Bond pointer to bond between \a *origin and \a *replacement
     165 * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
    166166 * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
    167167 * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
    … …  
    191191  InBondvector.SubtractVector(&TopOrigin->x);
    192192  bondlength = InBondvector.Norm();
    193    
     193
    194194   // is greater than typical bond distance? Then we have to correct periodically
    195195   // the problem is not the H being out of the box, but InBondvector have the wrong direction
    196    // due to TopReplacement or Origin being on the wrong side! 
    197   if (bondlength > BondDistance) { 
     196   // due to TopReplacement or Origin being on the wrong side!
     197  if (bondlength > BondDistance) {
    198198//    *out << Verbose(4) << "InBondvector is: ";
    199199//    InBondvector.Output(out);
    … …  
    215215//    *out << endl;
    216216  } // periodic correction finished
    217  
     217
    218218  InBondvector.Normalize();
    219219  // get typical bond length and store as scale factor for later
    … …  
    222222    cerr << Verbose(3) << "ERROR: There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
    223223    return false;
    224     BondRescale = bondlength; 
     224    BondRescale = bondlength;
    225225  } else {
    226226    if (!IsAngstroem)
    … …  
    273273      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
    274274//        *out << Verbose(3) << "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane." << endl;
    275        
     275
    276276        // determine the plane of these two with the *origin
    277277        AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
    … …  
    286286      Orthovector1.Normalize();
    287287      //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
    288      
     288
    289289      // create the two Hydrogens ...
    290290      FirstOtherAtom = new atom();
    … …  
    318318        SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle));
    319319      }
    320       FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance 
     320      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
    321321      SecondOtherAtom->x.Scale(&BondRescale);
    322322      //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
    323323      for(int i=NDIM;i--;) { // and make relative to origin atom
    324         FirstOtherAtom->x.x[i] += TopOrigin->x.x[i]; 
     324        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
    325325        SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
    326326      }
    … …  
    365365//      *out << endl;
    366366      AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
    367 //      *out << Verbose(3) << "Orthovector2: "; 
     367//      *out << Verbose(3) << "Orthovector2: ";
    368368//      Orthovector2.Output(out);
    369369//      *out << endl;
    370      
     370
    371371      // create correct coordination for the three atoms
    372372      alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
    … …  
    380380      factors[0] = d;
    381381      factors[1] = f;
    382       factors[2] = 0.; 
     382      factors[2] = 0.;
    383383      FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
    384384      factors[1] = -0.5*f;
    385       factors[2] = g; 
     385      factors[2] = g;
    386386      SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
    387       factors[2] = -g; 
     387      factors[2] = -g;
    388388      ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
    389389
    … …  
    437437 */
    438438bool molecule::AddXYZFile(string filename)
    439 { 
     439{
    440440  istringstream *input = NULL;
    441441  int NumberOfAtoms = 0; // atom number in xyz read
    … …  
    446446  string line;    // currently parsed line
    447447  double x[3];    // atom coordinates
    448  
     448
    449449  xyzfile.open(filename.c_str());
    450450  if (!xyzfile)
    … …  
    454454  input = new istringstream(line);
    455455  *input >> NumberOfAtoms;
    456   cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl; 
     456  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
    457457  getline(xyzfile,line,'\n'); // Read comment
    458458  cout << Verbose(1) << "Comment: " << line << endl;
    459  
     459
    460460  if (MDSteps == 0) // no atoms yet present
    461461    MDSteps++;
    … …  
    491491  xyzfile.close();
    492492  delete(input);
    493   return true; 
     493  return true;
    494494};
    495495
    … …  
    503503  atom *LeftAtom = NULL, *RightAtom = NULL;
    504504  atom *Walker = NULL;
    505  
     505
    506506  // copy all atoms
    507507  Walker = start;
    … …  
    510510    CurrentAtom = copy->AddCopyAtom(Walker);
    511511  }
    512  
     512
    513513  // copy all bonds
    514514  bond *Binder = first;
    … …  
    534534      copy->NoCyclicBonds++;
    535535    NewBond->Type = Binder->Type;
    536   } 
     536  }
    537537  // correct fathers
    538538  Walker = copy->start;
    … …  
    551551    copy->CreateListOfBondsPerAtom((ofstream *)&cout);
    552552  }
    553  
     553
    554554  return copy;
    555555};
    … …  
    576576
    577577/** Remove bond from bond chain list.
    578  * \todo Function not implemented yet 
     578 * \todo Function not implemented yet
    579579 * \param *pointer bond pointer
    580580 * \return true - bound found and removed, false - bond not found/removed
    … …  
    588588
    589589/** Remove every bond from bond chain list that atom \a *BondPartner is a constituent of.
    590  * \todo Function not implemented yet 
     590 * \todo Function not implemented yet
    591591 * \param *BondPartner atom to be removed
    592592 * \return true - bounds found and removed, false - bonds not found/removed
    … …  
    621621  Vector *min = new Vector;
    622622  Vector *max = new Vector;
    623  
     623
    624624  // gather min and max for each axis
    625625  ptr = start->next;  // start at first in list
    … …  
    667667{
    668668  Vector *min = new Vector;
    669  
     669
    670670//  *out << Verbose(3) << "Begin of CenterEdge." << endl;
    671671  atom *ptr = start->next;  // start at first in list
    … …  
    683683      }
    684684    }
    685 //    *out << Verbose(4) << "Maximum is "; 
     685//    *out << Verbose(4) << "Maximum is ";
    686686//    max->Output(out);
    687687//    *out << ", Minimum is ";
    … …  
    691691    max->AddVector(min);
    692692    Translate(min);
    693   } 
     693  }
    694694  delete(min);
    695695//  *out << Verbose(3) << "End of CenterEdge." << endl;
    696 }; 
     696};
    697697
    698698/** Centers the center of the atoms at (0,0,0).
    … …  
    704704  int Num = 0;
    705705  atom *ptr = start->next;  // start at first in list
    706  
     706
    707707  for(int i=NDIM;i--;) // zero center vector
    708708    center->x[i] = 0.;
    709    
     709
    710710  if (ptr != end) {   //list not empty?
    711711    while (ptr->next != end) {  // continue with second if present
    712712      ptr = ptr->next;
    713713      Num++;
    714       center->AddVector(&ptr->x);       
     714      center->AddVector(&ptr->x);
    715715    }
    716716    center->Scale(-1./Num); // divide through total number (and sign for direction)
    717717    Translate(center);
    718718  }
    719 }; 
     719};
    720720
    721721/** Returns vector pointing to center of gravity.
    … …  
    729729  Vector tmp;
    730730  double Num = 0;
    731  
     731
    732732  a->Zero();
    733733
    … …  
    737737      Num += 1.;
    738738      tmp.CopyVector(&ptr->x);
    739       a->AddVector(&tmp);       
     739      a->AddVector(&tmp);
    740740    }
    741741    a->Scale(-1./Num); // divide through total mass (and sign for direction)
    … …  
    757757        Vector tmp;
    758758  double Num = 0;
    759        
     759
    760760        a->Zero();
    761761
    … …  
    766766      tmp.CopyVector(&ptr->x);
    767767      tmp.Scale(ptr->type->mass);  // scale by mass
    768       a->AddVector(&tmp);       
     768      a->AddVector(&tmp);
    769769    }
    770770    a->Scale(-1./Num); // divide through total mass (and sign for direction)
    … …  
    789789    Translate(center);
    790790  }
    791 }; 
     791};
    792792
    793793/** Scales all atoms by \a *factor.
    … …  
    803803      Trajectories[ptr].R.at(j).Scale(factor);
    804804    ptr->x.Scale(factor);
    805   }     
    806 };
    807 
    808 /** Translate all atoms by given vector. 
     805  }
     806};
     807
     808/** Translate all atoms by given vector.
    809809 * \param trans[] translation vector.
    810810 */
    … …  
    818818      Trajectories[ptr].R.at(j).Translate(trans);
    819819    ptr->x.Translate(trans);
    820   }     
    821 };
    822 
    823 /** Mirrors all atoms against a given plane. 
     820  }
     821};
     822
     823/** Mirrors all atoms against a given plane.
    824824 * \param n[] normal vector of mirror plane.
    825825 */
    … …  
    833833      Trajectories[ptr].R.at(j).Mirror(n);
    834834    ptr->x.Mirror(n);
    835   }     
     835  }
    836836};
    837837
    … …  
    847847  bool flag;
    848848  Vector Testvector, Translationvector;
    849  
     849
    850850  do {
    851851    Center.Zero();
    … …  
    863863          if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    864864            for (int j=0;j<NDIM;j++) {
    865               tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j]; 
     865              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
    866866              if ((fabs(tmp)) > BondDistance) {
    867867                flag = false;
    … …  
    879879        cout << Verbose(1) << "vector is: ";
    880880        Testvector.Output((ofstream *)&cout);
    881         cout << endl;     
     881        cout << endl;
    882882#ifdef ADDHYDROGEN
    883883        // now also change all hydrogens
    … …  
    892892            cout << Verbose(1) << "Hydrogen vector is: ";
    893893            Testvector.Output((ofstream *)&cout);
    894             cout << endl;     
     894            cout << endl;
    895895          }
    896896        }
    … …  
    914914
    915915        CenterGravity(out, CenterOfGravity);
    916        
    917         // reset inertia tensor 
     916
     917        // reset inertia tensor
    918918        for(int i=0;i<NDIM*NDIM;i++)
    919919                InertiaTensor[i] = 0.;
    920        
     920
    921921        // sum up inertia tensor
    922922        while (ptr->next != end) {
    … …  
    943943        }
    944944        *out << endl;
    945        
     945
    946946        // diagonalize to determine principal axis system
    947947        gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
    … …  
    952952        gsl_eigen_symmv_free(T);
    953953        gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
    954        
     954
    955955        for(int i=0;i<NDIM;i++) {
    956956                *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
    957957                *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
    958958        }
    959        
     959
    960960        // check whether we rotate or not
    961961        if (DoRotate) {
    962           *out << Verbose(1) << "Transforming molecule into PAS ... "; 
     962          *out << Verbose(1) << "Transforming molecule into PAS ... ";
    963963          // the eigenvectors specify the transformation matrix
    964964          ptr = start;
    … …  
    972972
    973973          // summing anew for debugging (resulting matrix has to be diagonal!)
    974           // reset inertia tensor 
     974          // reset inertia tensor
    975975    for(int i=0;i<NDIM*NDIM;i++)
    976976      InertiaTensor[i] = 0.;
    977    
     977
    978978    // sum up inertia tensor
    979979    ptr = start;
    … …  
    10021002    *out << endl;
    10031003        }
    1004        
     1004
    10051005        // free everything
    10061006        delete(CenterOfGravity);
    … …  
    10311031
    10321032  CountElements();  // make sure ElementsInMolecule is up to date
    1033  
     1033
    10341034  // check file
    10351035  if (input == NULL) {
    … …  
    10871087              Trajectories[walker].U.at(MDSteps).x[d] += 0.5*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d]+Trajectories[walker].F.at(MDSteps).x[d])/IonMass;
    10881088            }
    1089 //            cout << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
     1089//            cout << "Integrated position&velocity of step " << (MDSteps) << ": (";
    10901090//            for (int d=0;d<NDIM;d++)
    10911091//              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
    … …  
    11211121  MDSteps++;
    11221122
    1123  
     1123
    11241124  // exit
    11251125  return true;
    11261126};
    11271127
    1128 /** Align all atoms in such a manner that given vector \a *n is along z axis. 
     1128/** Align all atoms in such a manner that given vector \a *n is along z axis.
    11291129 * \param n[] alignment vector.
    11301130 */
    … …  
    11451145    ptr = ptr->next;
    11461146    tmp = ptr->x.x[0];
    1147     ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2]; 
     1147    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
    11481148    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
    11491149    for (int j=0;j<MDSteps;j++) {
    11501150      tmp = Trajectories[ptr].R.at(j).x[0];
    1151       Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2]; 
     1151      Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
    11521152      Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
    11531153    }
    1154   }     
     1154  }
    11551155  // rotate n vector
    11561156  tmp = n->x[0];
    … …  
    11601160  n->Output((ofstream *)&cout);
    11611161  cout << endl;
    1162  
     1162
    11631163  // rotate on z-y plane
    11641164  ptr = start;
    … …  
    11681168    ptr = ptr->next;
    11691169    tmp = ptr->x.x[1];
    1170     ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2]; 
     1170    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
    11711171    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
    11721172    for (int j=0;j<MDSteps;j++) {
    11731173      tmp = Trajectories[ptr].R.at(j).x[1];
    1174       Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2]; 
     1174      Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
    11751175      Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
    11761176    }
    1177   }     
     1177  }
    11781178  // rotate n vector (for consistency check)
    11791179  tmp = n->x[1];
    11801180  n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
    11811181  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
    1182  
     1182
    11831183  cout << Verbose(1) << "alignment vector after second rotation: ";
    11841184  n->Output((ofstream *)&cout);
    … …  
    11911191 * \return true - succeeded, false - atom not found in list
    11921192 */
    1193 bool molecule::RemoveAtom(atom *pointer) 
    1194 { 
     1193bool molecule::RemoveAtom(atom *pointer)
     1194{
    11951195  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
    11961196    ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
    11971197  else
    1198     cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl; 
     1198    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
    11991199  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    12001200    ElementCount--;
    … …  
    12061206 * \return true - succeeded, false - atom not found in list
    12071207 */
    1208 bool molecule::CleanupMolecule() 
    1209 { 
    1210   return (cleanup(start,end) && cleanup(first,last)); 
     1208bool molecule::CleanupMolecule()
     1209{
     1210  return (cleanup(start,end) && cleanup(first,last));
    12111211};
    12121212
    … …  
    12221222  } else {
    12231223    cout << Verbose(0) << "Atom not found in list." << endl;
    1224     return NULL; 
     1224    return NULL;
    12251225  }
    12261226};
    … …  
    12711271  struct lsq_params *par = (struct lsq_params *)params;
    12721272  atom *ptr = par->mol->start;
    1273  
     1273
    12741274  // initialize vectors
    12751275  a.x[0] = gsl_vector_get(x,0);
    … …  
    13011301{
    13021302    int np = 6;
    1303    
     1303
    13041304   const gsl_multimin_fminimizer_type *T =
    13051305     gsl_multimin_fminimizer_nmsimplex;
    … …  
    13071307   gsl_vector *ss;
    13081308   gsl_multimin_function minex_func;
    1309  
     1309
    13101310   size_t iter = 0, i;
    13111311   int status;
    13121312   double size;
    1313  
     1313
    13141314   /* Initial vertex size vector */
    13151315   ss = gsl_vector_alloc (np);
    1316  
     1316
    13171317   /* Set all step sizes to 1 */
    13181318   gsl_vector_set_all (ss, 1.0);
    1319  
     1319
    13201320   /* Starting point */
    13211321   par->x = gsl_vector_alloc (np);
    13221322   par->mol = this;
    1323  
     1323
    13241324   gsl_vector_set (par->x, 0, 0.0);  // offset
    13251325   gsl_vector_set (par->x, 1, 0.0);
    … …  
    13281328   gsl_vector_set (par->x, 4, 0.0);
    13291329   gsl_vector_set (par->x, 5, 1.0);
    1330  
     1330
    13311331   /* Initialize method and iterate */
    13321332   minex_func.f = &LeastSquareDistance;
    13331333   minex_func.n = np;
    13341334   minex_func.params = (void *)par;
    1335  
     1335
    13361336   s = gsl_multimin_fminimizer_alloc (T, np);
    13371337   gsl_multimin_fminimizer_set (s, &minex_func, par->x, ss);
    1338  
     1338
    13391339   do
    13401340     {
    13411341       iter++;
    13421342       status = gsl_multimin_fminimizer_iterate(s);
    1343  
     1343
    13441344       if (status)
    13451345         break;
    1346  
     1346
    13471347       size = gsl_multimin_fminimizer_size (s);
    13481348       status = gsl_multimin_test_size (size, 1e-2);
    1349  
     1349
    13501350       if (status == GSL_SUCCESS)
    13511351         {
    13521352           printf ("converged to minimum at\n");
    13531353         }
    1354  
     1354
    13551355       printf ("%5d ", (int)iter);
    13561356       for (i = 0; i < (size_t)np; i++)
    … …  
    13611361     }
    13621362   while (status == GSL_CONTINUE && iter < 100);
    1363  
     1363
    13641364  for (i=0;i<(size_t)np;i++)
    13651365    gsl_vector_set(par->x, i, gsl_vector_get(s->x, i));
    … …  
    13781378  int ElementNo, AtomNo;
    13791379  CountElements();
    1380  
     1380
    13811381  if (out == NULL) {
    13821382    return false;
    … …  
    14131413  int ElementNo, AtomNo;
    14141414  CountElements();
    1415  
     1415
    14161416  if (out == NULL) {
    14171417    return false;
    … …  
    14601460  atom *Walker = start;
    14611461  while (Walker->next != end) {
    1462     Walker = Walker->next; 
     1462    Walker = Walker->next;
    14631463#ifdef ADDHYDROGEN
    14641464    if (Walker->type->Z != 1) {   // regard only non-hydrogen
    … …  
    14911491  int No = 0;
    14921492  time_t now;
    1493  
     1493
    14941494  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    14951495  walker = start;
    … …  
    15201520  int No = 0;
    15211521  time_t now;
    1522  
     1522
    15231523  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    15241524  walker = start;
    … …  
    15631563              Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
    15641564              if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
    1565                 NoNonHydrogen++; 
     1565                NoNonHydrogen++;
    15661566              Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
    15671567              Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
    … …  
    15711571            }
    15721572    } else
    1573         *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl; 
     1573        *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
    15741574  }
    15751575};
    … …  
    15831583        ElementsInMolecule[i] = 0;
    15841584        ElementCount = 0;
    1585        
     1585
    15861586  atom *walker = start;
    15871587  while (walker->next != end) {
    … …  
    16191619    Binder = Binder->next;
    16201620    if (Binder->Cyclic)
    1621       No++;   
     1621      No++;
    16221622  }
    16231623  delete(BackEdgeStack);
    … …  
    16791679 * Generally, we use the CSD approach to bond recognition, that is the the distance
    16801680 * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
    1681  * a threshold t = 0.4 Angstroem. 
     1681 * a threshold t = 0.4 Angstroem.
    16821682 * To make it O(N log N) the function uses the linked-cell technique as follows:
    16831683 * The procedure is step-wise:
    … …  
    16941694 * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
    16951695 */
     1696void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)// ofstream* out, double bonddistance, bool IsAngstroem)
     1697{
     1698
     1699        // 1 We will parse bonds out of the dbond file created by tremolo.
     1700                        int atom1, atom2, temp;
     1701                        atom *Walker, *OtherWalker;
     1702
     1703                if (!input)
     1704                {
     1705                        cout << Verbose(1) << "Opening silica failed \n";
     1706                };
     1707
     1708                        *input >> ws >> atom1;
     1709                        *input >> ws >> atom2;
     1710                cout << Verbose(1) << "Scanning file\n";
     1711                while (!input->eof()) // Check whether we read 2 items
     1712                {
     1713                                *input >> ws >> atom1;
     1714                                *input >> ws >> atom2;
     1715                        if(atom2<atom1) //Sort indizes of atoms in order
     1716                        {
     1717                                temp=atom1;
     1718                                atom1=atom2;
     1719                                atom2=temp;
     1720                        };
     1721
     1722                        Walker=start;
     1723                        while(Walker-> nr != atom1) // Find atom corresponding to first index
     1724                        {
     1725                                Walker = Walker->next;
     1726                        };
     1727                        OtherWalker = Walker->next;
     1728                        while(OtherWalker->nr != atom2) // Find atom corresponding to second index
     1729                        {
     1730                                OtherWalker= OtherWalker->next;
     1731                        };
     1732                        AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
     1733                        BondCount++; //Increase bond count of molecule
     1734                }
     1735
     1736                CreateListOfBondsPerAtom(out);
     1737
     1738};
    16961739void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
    16971740{
     1741
    16981742  atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
    16991743  int No, NoBonds, CandidateBondNo;
    … …  
    17041748  Vector x;
    17051749  int FalseBondDegree = 0;
    1706  
     1750
    17071751  BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
    17081752  *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
    … …  
    17111755    cleanup(first,last);
    17121756  }
    1713        
     1757
    17141758  // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
    17151759  CountAtoms(out);
    … …  
    17301774    for (int i=NumberCells;i--;)
    17311775      CellList[i] = NULL;
    1732  
     1776
    17331777    // 2b. put all atoms into its corresponding list
    17341778    Walker = start;
    … …  
    17511795      if (CellList[index] == NULL)  // allocate molecule if not done
    17521796        CellList[index] = new molecule(elemente);
    1753       OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference 
    1754       //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl; 
     1797      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
     1798      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
    17551799    }
    17561800    //for (int i=0;i<NumberCells;i++)
    17571801      //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
    1758      
     1802
     1803
    17591804    // 3a. go through every cell
    17601805    for (N[0]=divisor[0];N[0]--;)
    … …  
    17691814              Walker = Walker->next;
    17701815              //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
    1771               // 3c. check for possible bond between each atom in this and every one in the 27 cells 
     1816              // 3c. check for possible bond between each atom in this and every one in the 27 cells
    17721817              for (n[0]=-1;n[0]<=1;n[0]++)
    17731818                for (n[1]=-1;n[1]<=1;n[1]++)
    … …  
    18011846          }
    18021847        }
     1848
     1849
     1850
    18031851    // 4. free the cell again
    18041852    for (int i=NumberCells;i--;)
    … …  
    18071855      }
    18081856    Free((void **)&CellList, "molecule::CreateAdjacencyList - ** CellList");
    1809    
     1857
    18101858    // create the adjacency list per atom
    18111859    CreateListOfBondsPerAtom(out);
    1812                
     1860
    18131861    // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
    18141862    // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
    … …  
    18591907                        *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
    18601908          *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
    1861                
     1909
    18621910          // output bonds for debugging (if bond chain list was correctly installed)
    18631911          *out << Verbose(1) << endl << "From contents of bond chain list:";
    … …  
    18691917    *out << endl;
    18701918  } else
    1871         *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl; 
     1919        *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
    18721920  *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
    18731921  Free((void **)&matrix, "molecule::CreateAdjacencyList: *matrix");
    1874 };
     1922
     1923};
     1924
     1925
    18751926
    18761927/** Performs a Depth-First search on this molecule.
    … …  
    18931944  bond *Binder = NULL;
    18941945  bool BackStepping = false;
    1895  
     1946
    18961947  *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
    1897  
     1948
    18981949  ResetAllBondsToUnused();
    18991950  ResetAllAtomNumbers();
    … …  
    19081959    LeafWalker->Leaf = new molecule(elemente);
    19091960    LeafWalker->Leaf->AddCopyAtom(Root);
    1910    
     1961
    19111962    OldGraphNr = CurrentGraphNr;
    19121963    Walker = Root;
    … …  
    19191970          AtomStack->Push(Walker);
    19201971          CurrentGraphNr++;
    1921         }     
     1972        }
    19221973        do { // (3) if Walker has no unused egdes, go to (5)
    19231974          BackStepping = false; // reset backstepping flag for (8)
    … …  
    19532004          Binder = NULL;
    19542005      } while (1);  // (2)
    1955      
     2006
    19562007      // 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!
    19572008      if ((Walker == Root) && (Binder == NULL))
    19582009        break;
    1959        
    1960       // (5) if Ancestor of Walker is ... 
     2010
     2011      // (5) if Ancestor of Walker is ...
    19612012      *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
    19622013      if (Walker->Ancestor->GraphNr != Root->GraphNr) {
    … …  
    20012052        } while (OtherAtom != Walker);
    20022053        ComponentNumber++;
    2003    
     2054
    20042055        // (11) Root is separation vertex,  set Walker to Root and go to (4)
    20052056        Walker = Root;
    … …  
    20142065
    20152066    // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
    2016     *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;   
     2067    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
    20172068    LeafWalker->Leaf->Output(out);
    20182069    *out << endl;
    … …  
    20222073      //*out << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
    20232074      if (Root->GraphNr != -1) // if already discovered, step on
    2024         Root = Root->next; 
     2075        Root = Root->next;
    20252076    }
    20262077  }
    … …  
    20442095    *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
    20452096  }
    2046  
     2097
    20472098  *out << Verbose(1) << "Final graph info for each bond is:" << endl;
    20482099  Binder = first;
    … …  
    20552106    *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
    20562107    OutputComponentNumber(out, Binder->rightatom);
    2057     *out << ">." << endl; 
     2108    *out << ">." << endl;
    20582109    if (Binder->Cyclic) // cyclic ??
    20592110      *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
    … …  
    20702121 * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
    20712122 * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
    2072  * as cyclic and print out the cycles. 
     2123 * as cyclic and print out the cycles.
    20732124 * \param *out output stream for debugging
    20742125 * \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!
    … …  
    20812132  int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
    20822133  enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
    2083   class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring 
     2134  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring
    20842135  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
    20852136  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
    … …  
    20932144    ColorList[i] = white;
    20942145  }
    2095  
     2146
    20962147  *out << Verbose(1) << "Back edge list - ";
    20972148  BackEdgeStack->Output(out);
    2098  
     2149
    20992150  *out << Verbose(1) << "Analysing cycles ... " << endl;
    21002151  NumCycles = 0;
    … …  
    21022153    BackEdge = BackEdgeStack->PopFirst();
    21032154    // this is the target
    2104     Root = BackEdge->leftatom; 
     2155    Root = BackEdge->leftatom;
    21052156    // this is the source point
    2106     Walker = BackEdge->rightatom; 
     2157    Walker = BackEdge->rightatom;
    21072158    ShortestPathList[Walker->nr] = 0;
    21082159    BFSStack->ClearStack();  // start with empty BFS stack
    … …  
    21182169        if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
    21192170          OtherAtom = Binder->GetOtherAtom(Walker);
    2120 #ifdef ADDHYDROGEN         
     2171#ifdef ADDHYDROGEN
    21212172          if (OtherAtom->type->Z != 1) {
    21222173#endif
    … …  
    21272178              PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
    21282179              ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
    2129               *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; 
     2180              *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;
    21302181              //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
    21312182                *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
    … …  
    21372188            if (OtherAtom == Root)
    21382189              break;
    2139 #ifdef ADDHYDROGEN         
     2190#ifdef ADDHYDROGEN
    21402191          } else {
    21412192            *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
    … …  
    21752226      }
    21762227    } while ((!BFSStack->IsEmpty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
    2177    
     2228
    21782229    if (OtherAtom == Root) {
    21792230      // now climb back the predecessor list and thus find the cycle members
    … …  
    22032254      *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
    22042255    }
    2205    
     2256
    22062257    // now clean the lists
    22072258    while (!TouchedStack->IsEmpty()){
    … …  
    22132264  }
    22142265  if (MinRingSize != -1) {
    2215     // go over all atoms 
     2266    // go over all atoms
    22162267    Root = start;
    22172268    while(Root->next != end) {
    22182269      Root = Root->next;
    2219      
     2270
    22202271      if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
    22212272        Walker = Root;
    … …  
    22542305          }
    22552306          ColorList[Walker->nr] = black;
    2256           //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
     2307          //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
    22572308        }
    2258    
     2309
    22592310        // now clean the lists
    22602311        while (!TouchedStack->IsEmpty()){
    … …  
    23052356void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
    23062357{
    2307   for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++) 
     2358  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
    23082359    *out << vertex->ComponentNr[i] << "  ";
    23092360};
    … …  
    23832434{
    23842435  int c = 0;
    2385   int FragmentCount; 
     2436  int FragmentCount;
    23862437  // get maximum bond degree
    23872438  atom *Walker = start;
    … …  
    23932444  *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
    23942445  return FragmentCount;
    2395 }; 
     2446};
    23962447
    23972448/** Scans a single line for number and puts them into \a KeySet.
    23982449 * \param *out output stream for debugging
    23992450 * \param *buffer buffer to scan
    2400  * \param &CurrentSet filled KeySet on return 
     2451 * \param &CurrentSet filled KeySet on return
    24012452 * \return true - at least one valid atom id parsed, false - CurrentSet is empty
    24022453 */
    … …  
    24062457  int AtomNr;
    24072458  int status = 0;
    2408  
     2459
    24092460  line.str(buffer);
    24102461  while (!line.eof()) {
    … …  
    24422493  double TEFactor;
    24432494  char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
    2444  
     2495
    24452496  if (FragmentList == NULL) { // check list pointer
    24462497    FragmentList = new Graph;
    24472498  }
    2448  
     2499
    24492500  // 1st pass: open file and read
    24502501  *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
    … …  
    24752526    status = false;
    24762527  }
    2477  
     2528
    24782529  // 2nd pass: open TEFactors file and read
    24792530  *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
    … …  
    24872538        InputFile >> TEFactor;
    24882539        (*runner).second.second = TEFactor;
    2489         *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl; 
     2540        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
    24902541      } else {
    24912542        status = false;
    … …  
    25282579  if(output != NULL) {
    25292580    for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
    2530       for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) { 
     2581      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    25312582        if (sprinter != (*runner).first.begin())
    25322583          output << "\t";
    … …  
    25942645    status = false;
    25952646  }
    2596  
     2647
    25972648  return status;
    25982649};
    … …  
    26032654 * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
    26042655 * \return true - structure is equal, false - not equivalence
    2605  */ 
     2656 */
    26062657bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
    26072658{
    … …  
    26102661  bool status = true;
    26112662  char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
    2612  
     2663
    26132664  filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
    26142665  File.open(filename.str().c_str(), ios::out);
    … …  
    26692720  *out << endl;
    26702721  Free((void **)&buffer, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
    2671  
     2722
    26722723  return status;
    26732724};
    … …  
    26912742  for(int i=AtomCount;i--;)
    26922743    AtomMask[i] = false;
    2693  
     2744
    26942745  if (Order < 0) { // adaptive increase of BondOrder per site
    26952746    if (AtomMask[AtomCount] == true)  // break after one step
    … …  
    27312782          line >> ws >> Value; // skip time entry
    27322783          line >> ws >> Value;
    2733           No -= 1;  // indices start at 1 in file, not 0 
     2784          No -= 1;  // indices start at 1 in file, not 0
    27342785          //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
    27352786
    … …  
    27402791            // 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
    27412792            pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
    2742             map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first; 
     2793            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
    27432794            if (!InsertedElement.second) { // this root is already present
    2744               if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term 
    2745                 //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) 
     2795              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
     2796                //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)
    27462797                {  // if value is smaller, update value and order
    27472798                (*PresentItem).second.first = fabs(Value);
    … …  
    27812832        Walker = FindAtom(No);
    27822833        //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
    2783           *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl; 
     2834          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
    27842835          AtomMask[No] = true;
    27852836          status = true;
    27862837        //} else
    2787           //*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; 
     2838          //*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;
    27882839      }
    27892840      // close and done
    … …  
    28192870    if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
    28202871      status = true;
    2821      
     2872
    28222873    if (!status) {
    28232874      if (Order == 0)
    … …  
    28272878    }
    28282879  }
    2829  
     2880
    28302881  // print atom mask for debugging
    28312882  *out << "              ";
    … …  
    28362887    *out << (AtomMask[i] ? "t" : "f");
    28372888  *out << endl;
    2838  
     2889
    28392890  return status;
    28402891};
    … …  
    28502901  int AtomNo = 0;
    28512902  atom *Walker = NULL;
    2852  
     2903
    28532904  if (SortIndex != NULL) {
    28542905    *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
    … …  
    29082959  atom **ListOfAtoms = NULL;
    29092960  atom ***ListOfLocalAtoms = NULL;
    2910   bool *AtomMask = NULL; 
    2911  
     2961  bool *AtomMask = NULL;
     2962
    29122963  *out << endl;
    29132964#ifdef ADDHYDROGEN
    … …  
    29182969
    29192970  // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
    2920  
     2971
    29212972  // ===== 1. Check whether bond structure is same as stored in files ====
    2922  
     2973
    29232974  // fill the adjacency list
    29242975  CreateListOfBondsPerAtom(out);
    … …  
    29262977  // create lookup table for Atom::nr
    29272978  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);
    2928  
     2979
    29292980  // === compare it with adjacency file ===
    2930   FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms); 
     2981  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
    29312982  Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
    29322983
    … …  
    29573008    delete(LocalBackEdgeStack);
    29583009  }
    2959  
     3010
    29603011  // ===== 3. if structure still valid, parse key set file and others =====
    29613012  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
    … …  
    29633014  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
    29643015  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
    2965  
    2966   // =================================== Begin of FRAGMENTATION =============================== 
    2967   // ===== 6a. assign each keyset to its respective subgraph ===== 
     3016
     3017  // =================================== Begin of FRAGMENTATION ===============================
     3018  // ===== 6a. assign each keyset to its respective subgraph =====
    29683019  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
    29693020
    … …  
    29763027    FragmentationToDo = FragmentationToDo || CheckOrder;
    29773028    AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    2978     // ===== 6b. fill RootStack for each subgraph (second adaptivity check) ===== 
     3029    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    29793030    Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
    29803031
    … …  
    30003051  delete(ParsedFragmentList);
    30013052  delete[](MinimumRingSize);
    3002  
     3053
    30033054
    30043055  // ==================================== End of FRAGMENTATION ============================================
    … …  
    30063057  // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
    30073058  Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
    3008  
     3059
    30093060  // free subgraph memory again
    30103061  FragmentCounter = 0;
    … …  
    30313082    }
    30323083    *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
    3033    
     3084
    30343085    // ===== 9. Save fragments' configuration and keyset files et al to disk ===
    30353086    if (BondFragments->NumberOfMolecules != 0) {
    30363087      // create the SortIndex from BFS labels to order in the config file
    30373088      CreateMappingLabelsToConfigSequence(out, SortIndex);
    3038      
     3089
    30393090      *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
    30403091      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
    … …  
    30423093      else
    30433094        *out << Verbose(1) << "Some config writing failed." << endl;
    3044  
     3095
    30453096      // store force index reference file
    30463097      BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
    3047      
    3048       // store keysets file 
     3098
     3099      // store keysets file
    30493100      StoreKeySetFile(out, TotalGraph, configuration->configpath);
    3050  
    3051       // store Adjacency file 
     3101
     3102      // store Adjacency file
    30523103      StoreAdjacencyToFile(out, configuration->configpath);
    3053  
     3104
    30543105      // store Hydrogen saturation correction file
    30553106      BondFragments->AddHydrogenCorrection(out, configuration->configpath);
    3056      
     3107
    30573108      // store adaptive orders into file
    30583109      StoreOrderAtSiteFile(out, configuration->configpath);
    3059      
     3110
    30603111      // restore orbital and Stop values
    30613112      CalculateOrbitals(*configuration);
    3062      
     3113
    30633114      // free memory for bond part
    30643115      *out << Verbose(1) << "Freeing bond memory" << endl;
    30653116      delete(FragmentList); // remove bond molecule from memory
    3066       Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex"); 
     3117      Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
    30673118    } else
    30683119      *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
    3069   //} else 
     3120  //} else
    30703121  //  *out << Verbose(1) << "No fragments to store." << endl;
    30713122  *out << Verbose(0) << "End of bond fragmentation." << endl;
    … …  
    30933144  atom *Walker = NULL, *OtherAtom = NULL;
    30943145  ReferenceStack->Push(Binder);
    3095  
     3146
    30963147  do {  // go through all bonds and push local ones
    30973148    Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
    3098     if (Walker != NULL) // if this Walker exists in the subgraph ... 
     3149    if (Walker != NULL) // if this Walker exists in the subgraph ...
    30993150        for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
    31003151        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
    … …  
    31093160    ReferenceStack->Push(Binder);
    31103161  } while (FirstBond != Binder);
    3111  
     3162
    31123163  return status;
    31133164};
    … …  
    31853236      Walker->AdaptiveOrder = OrderArray[Walker->nr];
    31863237      Walker->MaxOrder = MaxArray[Walker->nr];
    3187       *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl; 
     3238      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl;
    31883239    }
    31893240    file.close();
    … …  
    31963247  Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
    31973248  Free((void **)&MaxArray, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
    3198  
     3249
    31993250  *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
    32003251  return status;
    … …  
    32543305  Walker = start;
    32553306  while (Walker->next != end) {
    3256     Walker = Walker->next; 
     3307    Walker = Walker->next;
    32573308    *out << Verbose(4) << "Atom " << Walker->Name << "/" << Walker->nr << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
    32583309    TotalDegree = 0;
    … …  
    32623313    }
    32633314    *out << " -- TotalDegree: " << TotalDegree << endl;
    3264   }     
     3315  }
    32653316  *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
    32663317};
    … …  
    32683319/** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
    32693320 * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
    3270  * white and putting into queue. 
     3321 * white and putting into queue.
    32713322 * \param *out output stream for debugging
    32723323 * \param *Mol Molecule class to add atoms to
    … …  
    32773328 * \param BondOrder maximum distance for vertices to add
    32783329 * \param IsAngstroem lengths are in angstroem or bohrradii
    3279  */ 
     3330 */
    32803331void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
    32813332{
    … …  
    33033354  }
    33043355  ShortestPathList[Root->nr] = 0;
    3305  
     3356
    33063357  // and go on ... Queue always contains all lightgray vertices
    33073358  while (!AtomStack->IsEmpty()) {
    … …  
    33113362    // followed by n+1 till top of stack.
    33123363    Walker = AtomStack->PopFirst(); // pop oldest added
    3313     *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl; 
     3364    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
    33143365    for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
    33153366      Binder = ListOfBondsPerAtom[Walker->nr][i];
    … …  
    33183369        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
    33193370        if (ColorList[OtherAtom->nr] == white) {
    3320           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) 
     3371          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)
    33213372            ColorList[OtherAtom->nr] = lightgray;
    33223373          PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
    33233374          ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
    3324           *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; 
     3375          *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;
    33253376          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
    33263377            *out << Verbose(3);
    … …  
    33703421          // This has to be a cyclic bond, check whether it's present ...
    33713422          if (AddedBondList[Binder->nr] == NULL) {
    3372             if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) { 
     3423            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
    33733424              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
    33743425              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
    … …  
    33853436    }
    33863437    ColorList[Walker->nr] = black;
    3387     *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
     3438    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
    33883439  }
    33893440  Free((void **)&PredecessorList, "molecule::BreadthFirstSearchAdd: **PredecessorList");
    … …  
    34093460
    34103461  *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
    3411  
     3462
    34123463  // reset parent list
    34133464  *out << Verbose(3) << "Resetting ParentList." << endl;
    34143465  for (int i=Father->AtomCount;i--;)
    34153466    ParentList[i] = NULL;
    3416  
     3467
    34173468  // fill parent list with sons
    34183469  *out << Verbose(3) << "Filling Parent List." << endl;
    … …  
    34553506 * \param *&Leaf KeySet to look through
    34563507 * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
    3457  * \param index of the atom suggested for removal 
     3508 * \param index of the atom suggested for removal
    34583509 */
    34593510int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
    … …  
    34613512  atom *Runner = NULL;
    34623513  int SP, Removal;
    3463  
     3514
    34643515  *out << Verbose(2) << "Looking for removal candidate." << endl;
    34653516  SP = -1; //0;  // not -1, so that Root is never removed
    … …  
    34793530/** Stores a fragment from \a KeySet into \a molecule.
    34803531 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
    3481  * molecule and adds missing hydrogen where bonds were cut. 
     3532 * molecule and adds missing hydrogen where bonds were cut.
    34823533 * \param *out output stream for debugging messages
    3483  * \param &Leaflet pointer to KeySet structure 
     3534 * \param &Leaflet pointer to KeySet structure
    34843535 * \param IsAngstroem whether we have Ansgtroem or bohrradius
    34853536 * \return pointer to constructed molecule
    … …  
    34923543  bool LonelyFlag = false;
    34933544  int size;
    3494  
     3545
    34953546//  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
    3496  
     3547
    34973548  Leaf->BondDistance = BondDistance;
    34983549  for(int i=NDIM*2;i--;)
    3499     Leaf->cell_size[i] = cell_size[i]; 
     3550    Leaf->cell_size[i] = cell_size[i];
    35003551
    35013552  // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
    … …  
    35103561    size++;
    35113562  }
    3512  
     3563
    35133564  // create the bonds between all: Make it an induced subgraph and add hydrogen
    35143565//  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
    … …  
    35203571    if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
    35213572      // create all bonds
    3522       for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father 
     3573      for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
    35233574        OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
    35243575//        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
    35253576        if (SonList[OtherFather->nr] != NULL) {
    3526 //          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl; 
     3577//          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
    35273578          if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
    35283579//            *out << Verbose(3) << "Adding Bond: ";
    3529 //            *out << 
     3580//            *out <<
    35303581            Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
    35313582//            *out << "." << endl;
    35323583            //NumBonds[Runner->nr]++;
    3533           } else { 
     3584          } else {
    35343585//            *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
    35353586          }
    35363587          LonelyFlag = false;
    35373588        } else {
    3538 //          *out << ", who has no son in this fragment molecule." << endl; 
     3589//          *out << ", who has no son in this fragment molecule." << endl;
    35393590#ifdef ADDHYDROGEN
    35403591          //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
    … …  
    35543605    while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
    35553606      Runner = Runner->next;
    3556 #endif       
     3607#endif
    35573608  }
    35583609  Leaf->CreateListOfBondsPerAtom(out);
    … …  
    35873638  StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
    35883639  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!
    3589   MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL; 
     3640  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
    35903641  MoleculeListClass *FragmentList = NULL;
    35913642  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
    … …  
    36373688                // clear snake stack
    36383689          SnakeStack->ClearStack();
    3639     //SnakeStack->TestImplementation(out, start->next);   
     3690    //SnakeStack->TestImplementation(out, start->next);
    36403691
    36413692    ///// INNER LOOP ////////////
    … …  
    36583709        }
    36593710        if (ShortestPathList[Walker->nr] == -1) {
    3660           ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1; 
     3711          ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1;
    36613712          TouchedStack->Push(Walker); // mark every atom for lists cleanup later, whose shortest path has been changed
    36623713        }
    … …  
    36963747                                OtherAtom = Binder->GetOtherAtom(Walker);
    36973748            if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
    3698               *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl; 
     3749              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
    36993750              //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
    37003751            } else { // otherwise check its colour and element
    37013752                                if (
    37023753#ifdef ADDHYDROGEN
    3703               (OtherAtom->type->Z != 1) && 
     3754              (OtherAtom->type->Z != 1) &&
    37043755#endif
    37053756                    (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
    … …  
    37113762                //} else {
    37123763                //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
    3713                 //} 
     3764                //}
    37143765                Walker = OtherAtom;
    37153766                break;
    37163767              } else {
    3717                 if (OtherAtom->type->Z == 1) 
     3768                if (OtherAtom->type->Z == 1)
    37183769                  *out << "Links to a hydrogen atom." << endl;
    3719                 else                 
     3770                else
    37203771                  *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
    37213772              }
    … …  
    37273778          *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
    37283779        }
    3729                         if (Walker != OtherAtom) {      // if no white neighbours anymore, color it black 
     3780                        if (Walker != OtherAtom) {      // if no white neighbours anymore, color it black
    37303781          *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
    37313782                                ColorVertexList[Walker->nr] = black;
    … …  
    37343785      }
    37353786        } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
    3736     *out << Verbose(2) << "Inner Looping is finished." << endl;   
     3787    *out << Verbose(2) << "Inner Looping is finished." << endl;
    37373788
    37383789    // if we reset all AtomCount atoms, we have again technically O(N^2) ...
    … …  
    37503801    }
    37513802  }
    3752   *out << Verbose(1) << "Outer Looping over all vertices is done." << endl; 
     3803  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
    37533804
    37543805  // copy together
    3755   *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl; 
     3806  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
    37563807  FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
    37573808  RunningIndex = 0;
    … …  
    38243875
    38253876  NumCombinations = 1 << SetDimension;
    3826  
     3877
    38273878  // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
    38283879  // von Endstuecken (aus den Bonds) hinzugefￜￜgt werden und fￜￜr verbleibende ANOVAOrder
    38293880  // rekursiv GraphCrawler in der nￜￜchsten Ebene aufgerufen werden
    3830  
     3881
    38313882  *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
    38323883  *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;
    38333884
    3834   // initialised touched list (stores added atoms on this level) 
     3885  // initialised touched list (stores added atoms on this level)
    38353886  *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
    38363887  for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
    38373888    TouchedList[TouchedIndex] = -1;
    38383889  TouchedIndex = 0;
    3839  
     3890
    38403891  // create every possible combination of the endpieces
    38413892  *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
    … …  
    38453896    for (int j=SetDimension;j--;)
    38463897      bits += (i & (1 << j)) >> j;
    3847      
     3898
    38483899    *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
    38493900    if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
    … …  
    38533904        bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
    38543905        if (bit) {  // if bit is set, we add this bond partner
    3855                 OtherWalker = BondsSet[j]->rightatom;   // rightatom is always the one more distant, i.e. the one to add 
     3906                OtherWalker = BondsSet[j]->rightatom;   // rightatom is always the one more distant, i.e. the one to add
    38563907          //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
    38573908          *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
    3858           TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr); 
     3909          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
    38593910          if (TestKeySetInsert.second) {
    38603911            TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
    … …  
    38693920        }
    38703921      }
    3871      
    3872       SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore 
     3922
     3923      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
    38733924      if (SpaceLeft > 0) {
    38743925        *out << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
    … …  
    38983949          }
    38993950          *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
    3900           SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits); 
     3951          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
    39013952          Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
    39023953        }
    … …  
    39073958        *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
    39083959        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
    3909           *out << (*runner) << " "; 
     3960          *out << (*runner) << " ";
    39103961        *out << endl;
    39113962        //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
    … …  
    39153966        //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
    39163967      }
    3917      
     3968
    39183969      // --3-- remove all added items in this level from snake stack
    39193970      *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
    … …  
    39263977      *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
    39273978      for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
    3928         *out << (*runner) << " "; 
     3979        *out << (*runner) << " ";
    39293980      *out << endl;
    39303981      TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
    … …  
    39333984    }
    39343985  }
    3935   Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList"); 
     3986  Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");
    39363987  *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
    39373988};
    … …  
    39423993 * \return true - connected, false - disconnected
    39433994 * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
    3944  */ 
     3995 */
    39453996bool molecule::CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment)
    39463997{
    … …  
    39483999  bool BondStatus = false;
    39494000  int size;
    3950  
     4001
    39514002  *out << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl;
    39524003  *out << Verbose(2) << "Disconnected atom: ";
    3953  
     4004
    39544005  // count number of atoms in graph
    39554006  size = 0;
    … …  
    39974048 * \param *out output stream for debugging
    39984049 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
    3999  * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on 
     4050 * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
    40004051 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
    40014052 * \return number of inserted fragments
    40024053 * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
    40034054 */
    4004 int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet) 
     4055int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
    40054056{
    40064057  int SP, AtomKeyNr;
    … …  
    40234074    FragmentSearch.BondsPerSPCount[i] = 0;
    40244075  FragmentSearch.BondsPerSPCount[0] = 1;
    4025   Binder = new bond(FragmentSearch.Root, FragmentSearch.Root); 
     4076  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
    40264077  add(Binder, FragmentSearch.BondsPerSPList[1]);
    4027  
     4078
    40284079  // do a BFS search to fill the SP lists and label the found vertices
    40294080  // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
    40304081  // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
    40314082  // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
    4032   // (EdgeinSPLevel) of this tree ... 
     4083  // (EdgeinSPLevel) of this tree ...
    40334084  // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
    40344085  // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
    … …  
    40834134    }
    40844135  }
    4085  
     4136
    40864137  // outputting all list for debugging
    40874138  *out << Verbose(0) << "Printing all found lists." << endl;
    … …  
    40924143      Binder = Binder->next;
    40934144      *out << Verbose(2) << *Binder << endl;
    4094     } 
    4095   }
    4096  
     4145    }
     4146  }
     4147
    40974148  // creating fragments with the found edge sets  (may be done in reverse order, faster)
    40984149  SP = -1;  // the Root <-> Root edge must be subtracted!
    … …  
    41014152    while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
    41024153      Binder = Binder->next;
    4103       SP ++; 
     4154      SP ++;
    41044155    }
    41054156  }
    … …  
    41084159    // start with root (push on fragment stack)
    41094160    *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
    4110     FragmentSearch.FragmentSet->clear(); 
     4161    FragmentSearch.FragmentSet->clear();
    41114162    *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
    41124163    // prepare the subset and call the generator
    41134164    BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
    41144165    BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond
    4115    
     4166
    41164167    SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
    4117    
     4168
    41184169    Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
    41194170  } else {
    … …  
    41244175  // remove root from stack
    41254176  *out << Verbose(0) << "Removing root again from stack." << endl;
    4126   FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);   
     4177  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
    41274178
    41284179  // free'ing the bonds lists
    … …  
    41434194  }
    41444195
    4145   // return list 
     4196  // return list
    41464197  *out << Verbose(0) << "End of PowerSetGenerator." << endl;
    41474198  return (FragmentSearch.FragmentCounter - Counter);
    … …  
    41744225    // remove bonds that are beyond bonddistance
    41754226    for(int i=NDIM;i--;)
    4176       Translationvector.x[i] = 0.; 
     4227      Translationvector.x[i] = 0.;
    41774228    // scan all bonds
    41784229    Binder = first;
    … …  
    42214272        }
    42224273      }
    4223       // re-add bond   
     4274      // re-add bond
    42244275      link(Binder, OtherBinder);
    42254276    } else {
    … …  
    42754326        IteratorB++;
    42764327      } // end of while loop
    4277     }// end of check in case of equal sizes 
     4328    }// end of check in case of equal sizes
    42784329  }
    42794330  return false; // if we reach this point, they are equal
    … …  
    43194370 * \param graph1 first (dest) graph
    43204371 * \param graph2 second (source) graph
    4321  * \param *counter keyset counter that gets increased 
     4372 * \param *counter keyset counter that gets increased
    43224373 */
    43234374inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
    … …  
    43644415  int RootKeyNr, RootNr;
    43654416  struct UniqueFragments FragmentSearch;
    4366  
     4417
    43674418  *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
    43684419
    … …  
    43874438    Walker = Walker->next;
    43884439    CompleteMolecule.insert(Walker->GetTrueFather()->nr);
    4389   } 
     4440  }
    43904441
    43914442  // 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
    … …  
    44014452    //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
    44024453    //  *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;
    4403     //} else 
     4454    //} else
    44044455    {
    44054456      // increase adaptive order by one
    44064457      Walker->GetTrueFather()->AdaptiveOrder++;
    44074458      Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
    4408  
     4459
    44094460      // initialise Order-dependent entries of UniqueFragments structure
    44104461      FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
    … …  
    44134464        FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
    44144465        FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
    4415         FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two 
     4466        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
    44164467        FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
    44174468        FragmentSearch.BondsPerSPCount[i] = 0;
    4418       } 
    4419  
     4469      }
     4470
    44204471      // allocate memory for all lower level orders in this 1D-array of ptrs
    44214472      NumLevels = 1 << (Order-1); // (int)pow(2,Order);
    … …  
    44234474      for (int i=0;i<NumLevels;i++)
    44244475        FragmentLowerOrdersList[RootNr][i] = NULL;
    4425      
     4476
    44264477      // create top order where nothing is reduced
    44274478      *out << Verbose(0) << "==============================================================================================================" << endl;
    44284479      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl; // , NumLevels is " << NumLevels << "
    4429  
     4480
    44304481      // Create list of Graphs of current Bond Order (i.e. F_{ij})
    44314482      FragmentLowerOrdersList[RootNr][0] =  new Graph;
    … …  
    44404491        // we don't have to dive into suborders! These keysets are all already created on lower orders!
    44414492        // this was all ancient stuff, when we still depended on the TEFactors (and for those the suborders were needed)
    4442        
     4493
    44434494//        if ((NumLevels >> 1) > 0) {
    44444495//          // create lower order fragments
    … …  
    44474498//          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)
    44484499//            // step down to next order at (virtual) boundary of powers of 2 in array
    4449 //            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
     4500//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
    44504501//              Order--;
    44514502//            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
    … …  
    44544505//              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
    44554506//              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
    4456 //       
     4507//
    44574508//              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
    44584509//              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
    … …  
    44854536      RootStack.push_back(RootKeyNr); // put back on stack
    44864537      RootNr++;
    4487  
     4538
    44884539      // free Order-dependent entries of UniqueFragments structure for next loop cycle
    44894540      Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
    … …  
    44914542        delete(FragmentSearch.BondsPerSPList[2*i]);
    44924543        delete(FragmentSearch.BondsPerSPList[2*i+1]);
    4493       } 
     4544      }
    44944545      Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
    44954546    }
    … …  
    45024553  Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
    45034554  delete(FragmentSearch.FragmentSet);
    4504  
    4505   // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein) 
     4555
     4556  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
    45064557  // 5433222211111111
    45074558  // 43221111
    … …  
    45234574    RootKeyNr = RootStack.front();
    45244575    RootStack.pop_front();
    4525     Walker = FindAtom(RootKeyNr); 
     4576    Walker = FindAtom(RootKeyNr);
    45264577    NumLevels = 1 << (Walker->AdaptiveOrder - 1);
    45274578    for(int i=0;i<NumLevels;i++) {
    … …  
    45364587  Free((void **)&FragmentLowerOrdersList, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
    45374588  Free((void **)&NumMoleculesOfOrder, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
    4538  
     4589
    45394590  *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
    45404591};
    … …  
    45704621  atom *Walker = NULL;
    45714622  bool result = true; // status of comparison
    4572  
    4573   *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl; 
     4623
     4624  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
    45744625  /// first count both their atoms and elements and update lists thereby ...
    45754626  //*out << Verbose(0) << "Counting atoms, updating list" << endl;
    … …  
    46184669    if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
    46194670      *out << Verbose(4) << "Centers of gravity don't match." << endl;
    4620       result = false; 
    4621     }
    4622   }
    4623  
     4671      result = false;
     4672    }
     4673  }
     4674
    46244675  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
    46254676  if (result) {
    … …  
    46374688      OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
    46384689    }
    4639  
     4690
    46404691    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    46414692    *out << Verbose(5) << "Sorting distances" << endl;
    … …  
    46484699    for(int i=AtomCount;i--;)
    46494700      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
    4650    
     4701
    46514702    /// ... and compare them step by step, whether the difference is individiually(!) below \a threshold for all
    46524703    *out << Verbose(4) << "Comparing distances" << endl;
    … …  
    46594710    Free((void **)&PermMap, "molecule::IsEqualToWithinThreshold: *PermMap");
    46604711    Free((void **)&OtherPermMap, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
    4661      
     4712
    46624713    /// free memory
    46634714    Free((void **)&Distances, "molecule::IsEqualToWithinThreshold: Distances");
    … …  
    46674718      result = false;
    46684719    }
    4669   } 
     4720  }
    46704721  /// return pointer to map if all distances were below \a threshold
    46714722  *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
    … …  
    46764727    *out << Verbose(3) << "Result: Not equal." << endl;
    46774728    return NULL;
    4678   } 
     4729  }
    46794730};
    46804731
    … …  
    47314782 * \param *output output stream of temperature file
    47324783 * \return file written (true), failure on writing file (false)
    4733  */ 
     4784 */
    47344785bool molecule::OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output)
    47354786{
    … …  
    47394790  if (output == NULL)
    47404791    return false;
    4741   else 
     4792  else
    47424793    *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
    47434794  for (int step=startstep;step < endstep; step++) { // loop over all time steps
Note: See TracChangeset for help on using the changeset viewer.