Changes in src/molecule.cpp [49e1ae:a67d19]

File:

• src/molecule.cpp (modified) (33 diffs)

src/molecule.cpp

@@ -23 +23 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "World.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -45 +46 @@
   for(int i=MAX_ELEMENTS;i--;)
     ElementsInMolecule[i] = 0;
-  cell_size[0] = cell_size[2] = cell_size[5]= 20.;
-  cell_size[1] = cell_size[3] = cell_size[4]= 0.;
-  strcpy(name,"none");
+  strcpy(name,World::get()->DefaultName);
 };
 
@@ -159 +158 @@
   double *matrix = NULL;
   bond *Binder = NULL;
+  double * const cell_size = World::get()->cell_size;
 
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -194 +194 @@
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
-    eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
     return false;
     BondRescale = bondlength;
@@ -237 +237 @@
             SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else {
-            eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name;
+            DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name);
           }
         }
@@ -274 +274 @@
       bondangle = TopOrigin->type->HBondAngle[1];
       if (bondangle == -1) {
-        eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+        DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
         return false;
         bondangle = 0;
@@ -396 +396 @@
       break;
     default:
-      eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl;
+      DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl);
       AllWentWell = false;
       break;
@@ -429 +429 @@
   input = new istringstream(line);
   *input >> NumberOfAtoms;
-  Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
+  DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl);
   getline(xyzfile,line,'\n'); // Read comment
-  Log() << Verbose(1) << "Comment: " << line << endl;
+  DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl);
 
   if (MDSteps == 0) // no atoms yet present
@@ -447 +447 @@
     Walker->type = elemente->FindElement(shorthand);
     if (Walker->type == NULL) {
-      eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.";
+      DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
       Walker->type = elemente->FindElement(1);
     }
@@ -543 +543 @@
     add(Binder, last);
   } else {
-    eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl);
   }
   return Binder;
@@ -555 +555 @@
 bool molecule::RemoveBond(bond *pointer)
 {
-  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
   pointer->leftatom->RegisterBond(pointer);
   pointer->rightatom->RegisterBond(pointer);
@@ -569 +569 @@
 bool molecule::RemoveBonds(atom *BondPartner)
 {
-  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
   BondList::const_iterator ForeRunner;
   while (!BondPartner->ListOfBonds.empty()) {
@@ -603 +603 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
+  double * const cell_size = World::get()->cell_size;
   cell_size[0] = dim->x[0];
   cell_size[1] = 0.;
@@ -621 +622 @@
     AtomCount--;
   } else
-    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "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--;
@@ -639 +640 @@
     ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
   else
-    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "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--;
@@ -664 +665 @@
     return walker;
   } else {
-    Log() << Verbose(0) << "Atom not found in list." << endl;
+    DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
     return NULL;
   }
@@ -680 +681 @@
     //mol->Output((ofstream *)&cout);
     //Log() << Verbose(0) << "===============================================" << endl;
-    Log() << Verbose(0) << text;
+    DoLog(0) && (Log() << Verbose(0) << text);
     cin >> No;
     ion = this->FindAtom(No);
@@ -693 +694 @@
 bool molecule::CheckBounds(const Vector *x) const
 {
+  double * const cell_size = World::get()->cell_size;
   bool result = true;
   int j =-1;
@@ -768 +770 @@
 void molecule::OutputListOfBonds() const
 {
-  Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
+  DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
   ActOnAllAtoms (&atom::OutputBondOfAtom );
-  Log() << Verbose(0) << endl;
+  DoLog(0) && (Log() << Verbose(0) << endl);
 };
 
@@ -827 +829 @@
   }
   if ((AtomCount == 0) || (i != AtomCount)) {
-    Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
+    DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);
     AtomCount = i;
 
@@ -843 +845 @@
         Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
         sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
-        Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
+        DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl);
         i++;
       }
     } else
-      Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
+      DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);
   }
 };
@@ -907 +909 @@
   bool result = true; // status of comparison
 
-  Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
+  DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
   /// first count both their atoms and elements and update lists thereby ...
   //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
@@ -919 +921 @@
   if (result) {
     if (AtomCount != OtherMolecule->AtomCount) {
-      Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
       result = false;
     } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
@@ -926 +928 @@
   if (result) {
     if (ElementCount != OtherMolecule->ElementCount) {
-      Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
+      DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl);
       result = false;
     } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
@@ -938 +940 @@
     }
     if (flag < MAX_ELEMENTS) {
-      Log() << Verbose(4) << "ElementsInMolecule don't match." << endl;
+      DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl);
       result = false;
     } else Log() << Verbose(4) << "ElementsInMolecule match." << endl;
@@ -944 +946 @@
   /// then determine and compare center of gravity for each molecule ...
   if (result) {
-    Log() << Verbose(5) << "Calculating Centers of Gravity" << endl;
+    DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
     DeterminePeriodicCenter(CenterOfGravity);
     OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
-    Log() << Verbose(5) << "Center of Gravity: ";
+    DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: ");
     CenterOfGravity.Output();
-    Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";
+    DoLog(0) && (Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ");
     OtherCenterOfGravity.Output();
-    Log() << Verbose(0) << endl;
+    DoLog(0) && (Log() << Verbose(0) << endl);
     if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
-      Log() << Verbose(4) << "Centers of gravity don't match." << endl;
+      DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
       result = false;
     }
@@ -960 +962 @@
   /// ... then make a list with the euclidian distance to this center for each atom of both molecules
   if (result) {
-    Log() << Verbose(5) << "Calculating distances" << endl;
+    DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
     Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
     OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
@@ -967 +969 @@
 
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
-    Log() << Verbose(5) << "Sorting distances" << endl;
+    DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
     PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
     OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
@@ -973 +975 @@
     gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
     PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
-    Log() << Verbose(5) << "Combining Permutation Maps" << endl;
+    DoLog(5) && (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
-    Log() << Verbose(4) << "Comparing distances" << endl;
+    DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
     flag = 0;
     for (int i=0;i<AtomCount;i++) {
-      Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
+      DoLog(5) && (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;
@@ -997 +999 @@
   }
   /// return pointer to map if all distances were below \a threshold
-  Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
+  DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
   if (result) {
-    Log() << Verbose(3) << "Result: Equal." << endl;
+    DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
     return PermutationMap;
   } else {
-    Log() << Verbose(3) << "Result: Not equal." << endl;
+    DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
     return NULL;
   }
@@ -1017 +1019 @@
 {
   atom *Walker = NULL, *OtherWalker = NULL;
-  Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
+  DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
   int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
   for (int i=AtomCount;i--;)
@@ -1024 +1026 @@
     for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
       AtomicMap[i] = i;
-    Log() << Verbose(4) << "Map is trivial." << endl;
+    DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
   } else {
-    Log() << Verbose(4) << "Map is ";
+    DoLog(4) && (Log() << Verbose(4) << "Map is ");
     Walker = start;
     while (Walker->next != end) {
@@ -1043 +1045 @@
         }
       }
-      Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";
-    }
-    Log() << Verbose(0) << endl;
-  }
-  Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl;
+      DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");
+    }
+    DoLog(0) && (Log() << Verbose(0) << endl);
+  }
+  DoLog(3) && (Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl);
   return AtomicMap;
 };