Changes in src/molecule.cpp [ead4e6:a67d19]

File:

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

src/molecule.cpp

@@ -6 +6 @@
 
 #include <cstring>
-#include <boost/bind.hpp>
-
-#include "World.hpp"
+
 #include "atom.hpp"
 #include "bond.hpp"
@@ -25 +23 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "World.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -31 +30 @@
  * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
  */
-molecule::molecule(const periodentafel * const teil) : elemente(teil), start(World::getInstance().createAtom()), end(World::getInstance().createAtom()),
+molecule::molecule(const periodentafel * const teil) : elemente(teil), start(new atom), end(new atom), 
   first(new bond(start, end, 1, -1)), last(new bond(start, end, 1, -1)), MDSteps(0), AtomCount(0), 
   BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0), NoCyclicBonds(0), BondDistance(0.), 
-  ActiveFlag(false), IndexNr(-1),
-  formula(this,boost::bind(&molecule::calcFormula,this)),
-  last_atom(0),
-  InternalPointer(start)
+  ActiveFlag(false), IndexNr(-1), last_atom(0), InternalPointer(start)
 {
   // init atom chain list
@@ -50 +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");
-};
-
-molecule *NewMolecule(){
-  return new molecule(World::getInstance().getPeriode());
-}
+  strcpy(name,World::get()->DefaultName);
+};
 
 /** Destructor of class molecule.
@@ -67 +57 @@
   delete(first);
   delete(last);
-  end->getWorld()->destroyAtom(end);
-  start->getWorld()->destroyAtom(start);
-};
-
-
-void DeleteMolecule(molecule *mol){
-  delete mol;
-}
-
-// getter and setter
-const std::string molecule::getName(){
-  return std::string(name);
-}
-
-void molecule::setName(const std::string _name){
-  OBSERVE;
-  strncpy(name,_name.c_str(),MAXSTRINGSIZE);
-}
-
-moleculeId_t molecule::getId(){
-  return id;
-}
-
-void molecule::setId(moleculeId_t _id){
-  id =_id;
-}
-
-const std::string molecule::getFormula(){
-  return *formula;
-}
-
-std::string molecule::calcFormula(){
-  std::map<atomicNumber_t,unsigned int> counts;
-  stringstream sstr;
-  periodentafel *periode = World::getInstance().getPeriode();
-  for(atom *Walker = start; Walker != end; Walker = Walker->next) {
-    counts[Walker->type->getNumber()]++;
-  }
-  std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
-  for(iter = counts.rbegin(); iter != counts.rend(); ++iter) {
-    atomicNumber_t Z = (*iter).first;
-    sstr << periode->FindElement(Z)->symbol << (*iter).second;
-  }
-  return sstr.str();
-}
+  delete(end);
+  delete(start);
+};
 
 
@@ -121 +69 @@
 bool molecule::AddAtom(atom *pointer)
 {
-  bool retval = false;
-  OBSERVE;
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
@@ -139 +85 @@
       }
     }
-    retval = add(pointer, end);
-  }
-  return retval;
+    return add(pointer, end);
+  } else
+    return false;
 };
 
@@ -151 +97 @@
 atom * molecule::AddCopyAtom(atom *pointer)
 {
-  atom *retval = NULL;
-  OBSERVE;
   if (pointer != NULL) {
-    atom *walker = pointer->clone();
+    atom *walker = new atom(pointer);
     walker->Name = Malloc<char>(strlen(pointer->Name) + 1, "atom::atom: *Name");
     strcpy (walker->Name, pointer->Name);
@@ -162 +106 @@
       NoNonHydrogen++;
     AtomCount++;
-    retval=walker;
-  }
-  return retval;
+    return walker;
+  } else
+    return NULL;
 };
 
@@ -203 +147 @@
 bool molecule::AddHydrogenReplacementAtom(bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
 {
-  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
-  OBSERVE;
   double bondlength;  // bond length of the bond to be replaced/cut
   double bondangle;  // bond angle of the bond to be replaced/cut
   double BondRescale;   // rescale value for the hydrogen bond length
+  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
   bond *FirstBond = NULL, *SecondBond = NULL; // Other bonds in double bond case to determine "other" plane
   atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
@@ -215 +158 @@
   double *matrix = NULL;
   bond *Binder = NULL;
+  double * const cell_size = World::get()->cell_size;
 
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -250 +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;
@@ -261 +205 @@
   switch(TopBond->BondDegree) {
     case 1:
-      FirstOtherAtom = World::getInstance().createAtom();    // new atom
+      FirstOtherAtom = new atom();    // new atom
       FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
       FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
@@ -293 +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);
           }
         }
@@ -318 +262 @@
 
       // create the two Hydrogens ...
-      FirstOtherAtom = World::getInstance().createAtom();
-      SecondOtherAtom = World::getInstance().createAtom();
+      FirstOtherAtom = new atom();
+      SecondOtherAtom = new atom();
       FirstOtherAtom->type = elemente->FindElement(1);
       SecondOtherAtom->type = elemente->FindElement(1);
@@ -330 +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;
@@ -373 +317 @@
     case 3:
       // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
-      FirstOtherAtom = World::getInstance().createAtom();
-      SecondOtherAtom = World::getInstance().createAtom();
-      ThirdOtherAtom = World::getInstance().createAtom();
+      FirstOtherAtom = new atom();
+      SecondOtherAtom = new atom();
+      ThirdOtherAtom = new atom();
       FirstOtherAtom->type = elemente->FindElement(1);
       SecondOtherAtom->type = elemente->FindElement(1);
@@ -452 +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;
@@ -469 +413 @@
 bool molecule::AddXYZFile(string filename)
 {
-
   istringstream *input = NULL;
   int NumberOfAtoms = 0; // atom number in xyz read
@@ -483 +426 @@
     return false;
 
-  OBSERVE;
   getline(xyzfile,line,'\n'); // Read numer of atoms in file
   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
     MDSteps++;
   for(i=0;i<NumberOfAtoms;i++){
-    Walker = World::getInstance().createAtom();
+    Walker = new atom;
     getline(xyzfile,line,'\n');
     istringstream *item = new istringstream(line);
@@ -505 +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);
     }
@@ -601 +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;
@@ -613 +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);
@@ -627 +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()) {
@@ -661 +603 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
+  double * const cell_size = World::get()->cell_size;
   cell_size[0] = dim->x[0];
   cell_size[1] = 0.;
@@ -679 +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--;
@@ -697 +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--;
@@ -722 +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;
   }
@@ -738 +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);
@@ -751 +694 @@
 bool molecule::CheckBounds(const Vector *x) const
 {
+  double * const cell_size = World::get()->cell_size;
   bool result = true;
   int j =-1;
@@ -826 +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);
 };
 
@@ -885 +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;
 
@@ -901 +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);
   }
 };
@@ -965 +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;
@@ -977 +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;
@@ -984 +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;
@@ -996 +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;
@@ -1002 +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;
     }
@@ -1018 +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");
@@ -1025 +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");
@@ -1031 +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;
@@ -1055 +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;
   }
@@ -1075 +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--;)
@@ -1082 +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) {
@@ -1101 +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;
 };
@@ -1141 +1085 @@
   }
 };
-
-void molecule::flipActiveFlag(){
-  ActiveFlag = !ActiveFlag;
-}