Changeset 543ce4 for molecuilder/src/molecule.cpp

Timestamp:

Nov 4, 2009, 7:56:04 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

4ef101, aa8542

Parents:

ec70ec

Message:

Huge change from ofstream * (const) out --> Log().

• first shift was done via regular expressions
• then via error messages from the code
• note that class atom, class element and class molecule kept in parts their output stream, was they print to file.
• make check runs fine
• MISSING: Verbosity is not fixed for everything (i.e. if no endl; is present and next has Verbose(0) ...)

Signed-off-by: Frederik Heber <heber@…>

File:

• molecuilder/src/molecule.cpp (modified) (53 diffs)

molecuilder/src/molecule.cpp

@@ -14 +14 @@
 #include "linkedcell.hpp"
 #include "lists.hpp"
+#include "log.hpp"
 #include "molecule.hpp"
 #include "memoryallocator.hpp"
@@ -143 +144 @@
  * \todo double and triple bonds splitting (always use the tetraeder angle!)
  */
-bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
+bool molecule::AddHydrogenReplacementAtom(bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
 {
   double bondlength;  // bond length of the bond to be replaced/cut
@@ -157 +158 @@
   bond *Binder = NULL;
 
-//  *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
+//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
   InBondvector.CopyVector(&TopReplacement->x);
@@ -167 +168 @@
    // due to TopReplacement or Origin being on the wrong side!
   if (bondlength > BondDistance) {
-//    *out << Verbose(4) << "InBondvector is: ";
+//    Log() << Verbose(4) << "InBondvector is: ";
 //    InBondvector.Output(out);
-//    *out << endl;
+//    Log() << Verbose(0) << endl;
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
@@ -182 +183 @@
     Free(&matrix);
     bondlength = InBondvector.Norm();
-//    *out << Verbose(4) << "Corrected InBondvector is now: ";
+//    Log() << Verbose(4) << "Corrected InBondvector is now: ";
 //    InBondvector.Output(out);
-//    *out << endl;
+//    Log() << Verbose(0) << endl;
   } // periodic correction finished
 
@@ -191 +192 @@
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
-    cerr << Verbose(3) << "ERROR: There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+    eLog() << Verbose(3) << "ERROR: There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
     return false;
     BondRescale = bondlength;
@@ -216 +217 @@
       FirstOtherAtom->x.AddVector(&InBondvector);  // ... and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
-//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
+//      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
 //      FirstOtherAtom->x.Output(out);
-//      *out << endl;
+//      Log() << Verbose(0) << endl;
       Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
       Binder->Cyclic = false;
@@ -234 +235 @@
             SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else {
-            *out << Verbose(3) << "WARNING: Detected more than four bonds for atom " << TopOrigin->Name;
+            Log() << Verbose(3) << "WARNING: Detected more than four bonds for atom " << TopOrigin->Name;
           }
         }
@@ -243 +244 @@
       }
       if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
-//        *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;
+//        Log() << 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;
 
         // determine the plane of these two with the *origin
@@ -250 +251 @@
         Orthovector1.GetOneNormalVector(&InBondvector);
       }
-      //*out << Verbose(3)<< "Orthovector1: ";
+      //Log() << Verbose(3)<< "Orthovector1: ";
       //Orthovector1.Output(out);
-      //*out << endl;
+      //Log() << Verbose(0) << endl;
       // orthogonal vector and bond vector between origin and replacement form the new plane
       Orthovector1.MakeNormalVector(&InBondvector);
       Orthovector1.Normalize();
-      //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
+      //Log() << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
 
       // create the two Hydrogens ...
@@ -271 +272 @@
       bondangle = TopOrigin->type->HBondAngle[1];
       if (bondangle == -1) {
-        *out << Verbose(3) << "ERROR: There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+        Log() << Verbose(3) << "ERROR: There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
         return false;
         bondangle = 0;
       }
       bondangle *= M_PI/180./2.;
-//      *out << Verbose(3) << "ReScaleCheck: InBondvector ";
+//      Log() << Verbose(3) << "ReScaleCheck: InBondvector ";
 //      InBondvector.Output(out);
-//      *out << endl;
-//      *out << Verbose(3) << "ReScaleCheck: Orthovector ";
+//      Log() << Verbose(0) << endl;
+//      Log() << Verbose(3) << "ReScaleCheck: Orthovector ";
 //      Orthovector1.Output(out);
-//      *out << endl;
-//      *out << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
+//      Log() << Verbose(0) << endl;
+//      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
       FirstOtherAtom->x.Zero();
       SecondOtherAtom->x.Zero();
@@ -291 +292 @@
       FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
       SecondOtherAtom->x.Scale(&BondRescale);
-      //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
+      //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
       for(int i=NDIM;i--;) { // and make relative to origin atom
         FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
@@ -299 +300 @@
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
       AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
-//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
+//      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
 //      FirstOtherAtom->x.Output(out);
-//      *out << endl;
-//      *out << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
+//      Log() << Verbose(0) << endl;
+//      Log() << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
 //      SecondOtherAtom->x.Output(out);
-//      *out << endl;
+//      Log() << Verbose(0) << endl;
       Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
       Binder->Cyclic = false;
@@ -332 +333 @@
       // we need to vectors orthonormal the InBondvector
       AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(&InBondvector);
-//      *out << Verbose(3) << "Orthovector1: ";
+//      Log() << Verbose(3) << "Orthovector1: ";
 //      Orthovector1.Output(out);
-//      *out << endl;
+//      Log() << Verbose(0) << endl;
       AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
-//      *out << Verbose(3) << "Orthovector2: ";
+//      Log() << Verbose(3) << "Orthovector2: ";
 //      Orthovector2.Output(out);
-//      *out << endl;
+//      Log() << Verbose(0) << endl;
 
       // create correct coordination for the three atoms
@@ -347 +348 @@
       f = b/sqrt(3.);   // length for Orthvector1
       g = b/2.;         // length for Orthvector2
-//      *out << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
-//      *out << Verbose(3) << "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << endl;
+//      Log() << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
+//      Log() << Verbose(3) << "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << endl;
       factors[0] = d;
       factors[1] = f;
@@ -373 +374 @@
       AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
       AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
-//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
+//      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
 //      FirstOtherAtom->x.Output(out);
-//      *out << endl;
-//      *out << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
+//      Log() << Verbose(0) << endl;
+//      Log() << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
 //      SecondOtherAtom->x.Output(out);
-//      *out << endl;
-//      *out << Verbose(4) << "Added " << *ThirdOtherAtom << " at: ";
+//      Log() << Verbose(0) << endl;
+//      Log() << Verbose(4) << "Added " << *ThirdOtherAtom << " at: ";
 //      ThirdOtherAtom->x.Output(out);
-//      *out << endl;
+//      Log() << Verbose(0) << endl;
       Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
       Binder->Cyclic = false;
@@ -393 +394 @@
       break;
     default:
-      cerr << "ERROR: BondDegree does not state single, double or triple bond!" << endl;
+      eLog() << Verbose(0) << "ERROR: BondDegree does not state single, double or triple bond!" << endl;
       AllWentWell = false;
       break;
   }
 
-//  *out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
+//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
   return AllWentWell;
 };
@@ -425 +426 @@
   input = new istringstream(line);
   *input >> NumberOfAtoms;
-  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
+  Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
   getline(xyzfile,line,'\n'); // Read comment
-  cout << Verbose(1) << "Comment: " << line << endl;
+  Log() << Verbose(1) << "Comment: " << line << endl;
 
   if (MDSteps == 0) // no atoms yet present
@@ -436 +437 @@
     istringstream *item = new istringstream(line);
     //istringstream input(line);
-    //cout << Verbose(1) << "Reading: " << line << endl;
+    //Log() << Verbose(1) << "Reading: " << line << endl;
     *item >> shorthand;
     *item >> x[0];
@@ -443 +444 @@
     Walker->type = elemente->FindElement(shorthand);
     if (Walker->type == NULL) {
-      cerr << "Could not parse the element at line: '" << line << "', setting to H.";
+      eLog() << Verbose(0) << "Could not parse the element at line: '" << line << "', setting to H.";
       Walker->type = elemente->FindElement(1);
     }
@@ -496 +497 @@
 
   // copy values
-  copy->CountAtoms((ofstream *)&cout);
+  copy->CountAtoms();
   copy->CountElements();
   if (first->next != last) {  // if adjaceny list is present
@@ -517 +518 @@
   ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
 
-  //TODO: copy->BuildInducedSubgraph((ofstream *)&cout, this);
+  //TODO: copy->BuildInducedSubgraph(this);
 
   return copy;
@@ -539 +540 @@
     add(Binder, last);
   } else {
-    cerr << Verbose(1) << "ERROR: Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
+    eLog() << Verbose(1) << "ERROR: Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
   }
   return Binder;
@@ -551 +552 @@
 bool molecule::RemoveBond(bond *pointer)
 {
-  //cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
   pointer->leftatom->RegisterBond(pointer);
   pointer->rightatom->RegisterBond(pointer);
@@ -565 +566 @@
 bool molecule::RemoveBonds(atom *BondPartner)
 {
-  //cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
   BondList::const_iterator ForeRunner;
   while (!BondPartner->ListOfBonds.empty()) {
@@ -617 +618 @@
     AtomCount--;
   } else
-    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+    eLog() << Verbose(0) << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
@@ -635 +636 @@
     ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
   else
-    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+    eLog() << Verbose(0) << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
@@ -657 +658 @@
   atom * walker = find(&Nr, start,end);
   if (walker != NULL) {
-    //cout << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
+    //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
     return walker;
   } else {
-    cout << Verbose(0) << "Atom not found in list." << endl;
+    Log() << Verbose(0) << "Atom not found in list." << endl;
     return NULL;
   }
@@ -673 +674 @@
   atom *ion = NULL;
   do {
-    //cout << Verbose(0) << "============Atom list==========================" << endl;
+    //Log() << Verbose(0) << "============Atom list==========================" << endl;
     //mol->Output((ofstream *)&cout);
-    //cout << Verbose(0) << "===============================================" << endl;
-    cout << Verbose(0) << text;
+    //Log() << Verbose(0) << "===============================================" << endl;
+    Log() << Verbose(0) << text;
     cin >> No;
     ion = this->FindAtom(No);
@@ -702 +703 @@
  * \param *out output stream
  */
-bool molecule::Output(ofstream *out)
+bool molecule::Output(ofstream * const output)
 {
   int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
@@ -711 +712 @@
     ElementNo[i] = 0;
   }
-  if (out == NULL) {
+  if (output == NULL) {
     return false;
   } else {
-    *out << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
+    *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
     SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
     int current=1;
@@ -721 +722 @@
         ElementNo[i] = current++;
     }
-    ActOnAllAtoms( &atom::OutputArrayIndexed, (ofstream *)out, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
+    ActOnAllAtoms( &atom::OutputArrayIndexed, output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
     return true;
   }
@@ -729 +730 @@
  * \param *out output stream
  */
-bool molecule::OutputTrajectories(ofstream *out)
+bool molecule::OutputTrajectories(ofstream * const output)
 {
   int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
   CountElements();
 
-  if (out == NULL) {
+  if (output == NULL) {
     return false;
   } else {
     for (int step = 0; step < MDSteps; step++) {
       if (step == 0) {
-        *out << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
+        *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
       } else {
-        *out << "# ====== MD step " << step << " =========" << endl;
+        *output << "# ====== MD step " << step << " =========" << endl;
       }
       for (int i=0;i<MAX_ELEMENTS;++i) {
@@ -753 +754 @@
           ElementNo[i] = current++;
       }
-      ActOnAllAtoms( &atom::OutputTrajectory, out, (const int *)ElementNo, AtomNo, (const int)step );
+      ActOnAllAtoms( &atom::OutputTrajectory, output, (const int *)ElementNo, AtomNo, (const int)step );
     }
     return true;
@@ -762 +763 @@
  * \param *out output stream
  */
-void molecule::OutputListOfBonds(ofstream *out) const
-{
-  *out << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
-  ActOnAllAtoms (&atom::OutputBondOfAtom, out);
-  *out << endl;
+void molecule::OutputListOfBonds() const
+{
+  Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
+  ActOnAllAtoms (&atom::OutputBondOfAtom );
+  Log() << Verbose(0) << endl;
 };
 
@@ -772 +773 @@
  * \param *out stream pointer
  */
-bool molecule::Checkout(ofstream *out)  const
-{
-  return elemente->Checkout(out, ElementsInMolecule);
+bool molecule::Checkout(ofstream * const output)  const
+{
+  return elemente->Checkout(output, ElementsInMolecule);
 };
 
@@ -780 +781 @@
  * \param *out output stream
  */
-bool molecule::OutputTrajectoriesXYZ(ofstream *out)
+bool molecule::OutputTrajectoriesXYZ(ofstream * const output)
 {
   time_t now;
 
-  if (out != NULL) {
+  if (output != NULL) {
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     for (int step=0;step<MDSteps;step++) {
-      *out << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
-      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, out, step );
+      *output << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
+      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     }
     return true;
@@ -798 +799 @@
 * \param *out output stream
  */
-bool molecule::OutputXYZ(ofstream *out) const
+bool molecule::OutputXYZ(ofstream * const output) const
 {
   time_t now;
 
-  if (out != NULL) {
+  if (output != NULL) {
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
-    *out << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
-    ActOnAllAtoms( &atom::OutputXYZLine, out );
+    *output << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
+    ActOnAllAtoms( &atom::OutputXYZLine, output );
     return true;
   } else
@@ -814 +815 @@
  * \param *out output stream for debugging
  */
-void molecule::CountAtoms(ofstream *out)
+void molecule::CountAtoms()
 {
   int i = 0;
@@ -823 +824 @@
   }
   if ((AtomCount == 0) || (i != AtomCount)) {
-    *out << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
+    Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
     AtomCount = i;
 
@@ -839 +840 @@
         Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
         sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
-        *out << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
+        Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
         i++;
       }
     } else
-      *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
+      Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
   }
 };
@@ -870 +871 @@
   for(int i=MAX_ELEMENTS;i--;) {
     if (ElementsInMolecule[i] != 0) {
-      //cout << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
+      //Log() << Verbose(0) << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
       configuration.MaxPsiDouble += ElementsInMolecule[i]*((int)elemente->FindElement(i)->Valence);
     }
@@ -894 +895 @@
  * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
  */
-int * molecule::IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold)
+int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)
 {
   int flag;
@@ -903 +904 @@
   bool result = true; // status of comparison
 
-  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
+  Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
   /// first count both their atoms and elements and update lists thereby ...
-  //*out << Verbose(0) << "Counting atoms, updating list" << endl;
-  CountAtoms(out);
-  OtherMolecule->CountAtoms(out);
+  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
+  CountAtoms();
+  OtherMolecule->CountAtoms();
   CountElements();
   OtherMolecule->CountElements();
@@ -915 +916 @@
   if (result) {
     if (AtomCount != OtherMolecule->AtomCount) {
-      *out << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+      Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
       result = false;
-    } else *out << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
   }
   /// -# ElementCount
   if (result) {
     if (ElementCount != OtherMolecule->ElementCount) {
-      *out << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
+      Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
       result = false;
-    } else *out << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
+    } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
   }
   /// -# ElementsInMolecule
   if (result) {
     for (flag=MAX_ELEMENTS;flag--;) {
-      //*out << Verbose(5) << "Element " <<  flag << ": " << ElementsInMolecule[flag] << " <-> " << OtherMolecule->ElementsInMolecule[flag] << "." << endl;
+      //Log() << Verbose(5) << "Element " <<  flag << ": " << ElementsInMolecule[flag] << " <-> " << OtherMolecule->ElementsInMolecule[flag] << "." << endl;
       if (ElementsInMolecule[flag] != OtherMolecule->ElementsInMolecule[flag])
         break;
     }
     if (flag < MAX_ELEMENTS) {
-      *out << Verbose(4) << "ElementsInMolecule don't match." << endl;
+      Log() << Verbose(4) << "ElementsInMolecule don't match." << endl;
       result = false;
-    } else *out << Verbose(4) << "ElementsInMolecule match." << endl;
+    } else Log() << Verbose(4) << "ElementsInMolecule match." << endl;
   }
   /// then determine and compare center of gravity for each molecule ...
   if (result) {
-    *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
+    Log() << Verbose(5) << "Calculating Centers of Gravity" << endl;
     DeterminePeriodicCenter(CenterOfGravity);
     OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
-    *out << Verbose(5) << "Center of Gravity: ";
-    CenterOfGravity.Output(out);
-    *out << endl << Verbose(5) << "Other Center of Gravity: ";
-    OtherCenterOfGravity.Output(out);
-    *out << endl;
+    Log() << Verbose(5) << "Center of Gravity: ";
+    CenterOfGravity.Output();
+    Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";
+    OtherCenterOfGravity.Output();
+    Log() << Verbose(0) << endl;
     if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
-      *out << Verbose(4) << "Centers of gravity don't match." << endl;
+      Log() << Verbose(4) << "Centers of gravity don't match." << endl;
       result = false;
     }
@@ -956 +957 @@
   /// ... then make a list with the euclidian distance to this center for each atom of both molecules
   if (result) {
-    *out << Verbose(5) << "Calculating distances" << endl;
+    Log() << Verbose(5) << "Calculating distances" << endl;
     Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
     OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
@@ -963 +964 @@
 
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
-    *out << Verbose(5) << "Sorting distances" << endl;
+    Log() << Verbose(5) << "Sorting distances" << endl;
     PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
     OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
@@ -969 +970 @@
     gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
     PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
-    *out << Verbose(5) << "Combining Permutation Maps" << endl;
+    Log() << Verbose(5) << "Combining Permutation Maps" << endl;
     for(int i=AtomCount;i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
 
     /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
-    *out << Verbose(4) << "Comparing distances" << endl;
+    Log() << Verbose(4) << "Comparing distances" << endl;
     flag = 0;
     for (int i=0;i<AtomCount;i++) {
-      *out << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
+      Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
       if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
         flag = 1;
@@ -993 +994 @@
   }
   /// return pointer to map if all distances were below \a threshold
-  *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
+  Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
   if (result) {
-    *out << Verbose(3) << "Result: Equal." << endl;
+    Log() << Verbose(3) << "Result: Equal." << endl;
     return PermutationMap;
   } else {
-    *out << Verbose(3) << "Result: Not equal." << endl;
+    Log() << Verbose(3) << "Result: Not equal." << endl;
     return NULL;
   }
@@ -1010 +1011 @@
  * \todo make this with a good sort O(n), not O(n^2)
  */
-int * molecule::GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule)
+int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
 {
   atom *Walker = NULL, *OtherWalker = NULL;
-  *out << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
+  Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
   int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
   for (int i=AtomCount;i--;)
@@ -1020 +1021 @@
     for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
       AtomicMap[i] = i;
-    *out << Verbose(4) << "Map is trivial." << endl;
+    Log() << Verbose(4) << "Map is trivial." << endl;
   } else {
-    *out << Verbose(4) << "Map is ";
+    Log() << Verbose(4) << "Map is ";
     Walker = start;
     while (Walker->next != end) {
@@ -1034 +1035 @@
       //for (int i=0;i<AtomCount;i++) { // search atom
         //for (int j=0;j<OtherMolecule->AtomCount;j++) {
-          //*out << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
+          //Log() << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
           if (Walker->father == OtherWalker)
             AtomicMap[Walker->nr] = OtherWalker->nr;
         }
       }
-      *out << AtomicMap[Walker->nr] << "\t";
-    }
-    *out << endl;
-  }
-  *out << Verbose(3) << "End of GetFatherAtomicMap." << endl;
+      Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";
+    }
+    Log() << Verbose(0) << endl;
+  }
+  Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl;
   return AtomicMap;
 };
@@ -1050 +1051 @@
  * We simply use the formula equivaleting temperature and kinetic energy:
  * \f$k_B T = \sum_i m_i v_i^2\f$
- * \param *out output stream for debugging
+ * \param *output output stream of temperature file
  * \param startstep first MD step in molecule::Trajectories
  * \param endstep last plus one MD step in molecule::Trajectories
- * \param *output output stream of temperature file
  * \return file written (true), failure on writing file (false)
  */
-bool molecule::OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output)
+bool molecule::OutputTemperatureFromTrajectories(ofstream * const output, int startstep, int endstep)
 {
   double temperature;