Changes in src/molecule.cpp [bab12a:f8e486]

File:

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

src/molecule.cpp

@@ -4 +4 @@
  *
  */
+
+#include "Helpers/MemDebug.hpp"
 
 #include <cstring>
@@ -35 +37 @@
  * 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()),
-  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)
-{
-  // init atom chain list
-  start->father = NULL;
-  end->father = NULL;
-  link(start,end);
-
-  // init bond chain list
-  link(first,last);
+molecule::molecule(const periodentafel * const teil) :
+  Observable("molecule"),
+  elemente(teil),  MDSteps(0),  BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0),
+  NoCyclicBonds(0), BondDistance(0.),  ActiveFlag(false), IndexNr(-1),
+  formula(this,boost::bind(&molecule::calcFormula,this),"formula"),
+  AtomCount(this,boost::bind(&molecule::doCountAtoms,this),"AtomCount"), last_atom(0),  InternalPointer(begin())
+{
 
   // other stuff
@@ -67 +61 @@
 {
   CleanupMolecule();
-  delete(first);
-  delete(last);
-  end->getWorld()->destroyAtom(end);
-  start->getWorld()->destroyAtom(start);
 };
 
@@ -81 +71 @@
 const std::string molecule::getName(){
   return std::string(name);
+}
+
+int molecule::getAtomCount() const{
+  return *AtomCount;
 }
 
@@ -104 +98 @@
   stringstream sstr;
   periodentafel *periode = World::getInstance().getPeriode();
-  for(atom *Walker = start; Walker != end; Walker = Walker->next) {
-    counts[Walker->type->getNumber()]++;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    counts[(*iter)->type->getNumber()]++;
   }
   std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
@@ -115 +109 @@
 }
 
+/************************** Access to the List of Atoms ****************/
+
+
+molecule::iterator molecule::begin(){
+  return molecule::iterator(atoms.begin(),this);
+}
+
+molecule::const_iterator molecule::begin() const{
+  return atoms.begin();
+}
+
+molecule::iterator molecule::end(){
+  return molecule::iterator(atoms.end(),this);
+}
+
+molecule::const_iterator molecule::end() const{
+  return atoms.end();
+}
+
+bool molecule::empty() const
+{
+  return (begin() == end());
+}
+
+size_t molecule::size() const
+{
+  size_t counter = 0;
+  for (molecule::const_iterator iter = begin(); iter != end (); ++iter)
+    counter++;
+  return counter;
+}
+
+molecule::const_iterator molecule::erase( const_iterator loc )
+{
+  molecule::const_iterator iter = loc;
+  iter--;
+  atom* atom = *loc;
+  atoms.erase( loc );
+  atom->removeFromMolecule();
+  return iter;
+}
+
+molecule::const_iterator molecule::erase( atom * key )
+{
+  cout << "trying to erase atom" << endl;
+  molecule::const_iterator iter = find(key);
+  if (iter != end()){
+    atoms.erase( iter++ );
+    key->removeFromMolecule();
+  }
+  return iter;
+}
+
+molecule::const_iterator molecule::find ( atom * key ) const
+{
+  return atoms.find( key );
+}
+
+pair<molecule::iterator,bool> molecule::insert ( atom * const key )
+{
+  pair<atomSet::iterator,bool> res = atoms.insert(key);
+  return pair<iterator,bool>(iterator(res.first,this),res.second);
+}
+
+bool molecule::containsAtom(atom* key){
+  return atoms.count(key);
+}
 
 /** Adds given atom \a *pointer from molecule list.
@@ -123 +184 @@
 bool molecule::AddAtom(atom *pointer)
 {
-  bool retval = false;
   OBSERVE;
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    pointer->nr = last_atom++;  // increase number within molecule
-    AtomCount++;
     if (pointer->type != NULL) {
       if (ElementsInMolecule[pointer->type->Z] == 0)
@@ -141 +199 @@
       }
     }
-    retval = add(pointer, end);
-  }
-  return retval;
+    insert(pointer);
+    pointer->setMolecule(this);
+  }
+  return true;
 };
 
@@ -157 +216 @@
   if (pointer != NULL) {
     atom *walker = pointer->clone();
-    stringstream sstr;
-    sstr << pointer->getName();
-    walker->setName(sstr.str());
+    walker->setName(pointer->getName());
     walker->nr = last_atom++;  // increase number within molecule
-    add(walker, end);
+    insert(walker);
     if ((pointer->type != NULL) && (pointer->type->Z != 1))
       NoNonHydrogen++;
-    AtomCount++;
     retval=walker;
   }
@@ -242 +298 @@
     Orthovector1.MatrixMultiplication(matrix);
     InBondvector -= Orthovector1; // subtract just the additional translation
-    Free(&matrix);
+    delete[](matrix);
     bondlength = InBondvector.Norm();
 //    Log() << Verbose(4) << "Corrected InBondvector is now: ";
@@ -251 +307 @@
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
+  ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
@@ -472 +529 @@
       break;
   }
-  Free(&matrix);
+  delete[](matrix);
 
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
@@ -555 +612 @@
 
   // copy all bonds
-  bond *Binder = first;
+  bond *Binder = NULL;
   bond *NewBond = NULL;
-  while(Binder->next != last) {
-    Binder = Binder->next;
-
-    // get the pendant atoms of current bond in the copy molecule
-    copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->leftatom, (const atom **)&LeftAtom );
-    copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->rightatom, (const atom **)&RightAtom );
-
-    NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
-    NewBond->Cyclic = Binder->Cyclic;
-    if (Binder->Cyclic)
-      copy->NoCyclicBonds++;
-    NewBond->Type = Binder->Type;
-  }
+  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
+    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
+      if ((*BondRunner)->leftatom == *AtomRunner) {
+        Binder = (*BondRunner);
+
+        // get the pendant atoms of current bond in the copy molecule
+        copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->leftatom, (const atom **)&LeftAtom );
+        copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->rightatom, (const atom **)&RightAtom );
+
+        NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
+        NewBond->Cyclic = Binder->Cyclic;
+        if (Binder->Cyclic)
+          copy->NoCyclicBonds++;
+        NewBond->Type = Binder->Type;
+      }
   // correct fathers
   ActOnAllAtoms( &atom::CorrectFather );
 
   // copy values
-  copy->CountAtoms();
   copy->CountElements();
-  if (first->next != last) {  // if adjaceny list is present
+  if (hasBondStructure()) {  // if adjaceny list is present
     copy->BondDistance = BondDistance;
   }
@@ -608 +666 @@
 bond * molecule::AddBond(atom *atom1, atom *atom2, int degree)
 {
+  OBSERVE;
   bond *Binder = NULL;
-  if ((atom1 != NULL) && (FindAtom(atom1->nr) != NULL) && (atom2 != NULL) && (FindAtom(atom2->nr) != NULL)) {
-    Binder = new bond(atom1, atom2, degree, BondCount++);
-    atom1->RegisterBond(Binder);
-    atom2->RegisterBond(Binder);
-    if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
-      NoNonBonds++;
-    add(Binder, last);
-  } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Could not add bond between " << atom1->getName() << " and " << atom2->getName() << " as one or both are not present in the molecule." << endl);
-  }
+
+  // some checks to make sure we are able to create the bond
+  ASSERT(atom1, "First atom in bond-creation was an invalid pointer");
+  ASSERT(atom2, "Second atom in bond-creation was an invalid pointer");
+  ASSERT(FindAtom(atom1->nr),"First atom in bond-creation was not part of molecule");
+  ASSERT(FindAtom(atom2->nr),"Second atom in bond-creation was not parto of molecule");
+
+  Binder = new bond(atom1, atom2, degree, BondCount++);
+  atom1->RegisterBond(Binder);
+  atom2->RegisterBond(Binder);
+  if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
+    NoNonBonds++;
+
   return Binder;
 };
@@ -630 +692 @@
 {
   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
-  pointer->leftatom->RegisterBond(pointer);
-  pointer->rightatom->RegisterBond(pointer);
-  removewithoutcheck(pointer);
+  delete(pointer);
   return true;
 };
@@ -692 +752 @@
 bool molecule::RemoveAtom(atom *pointer)
 {
+  ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
+  OBSERVE;
   if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
     ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
-    AtomCount--;
   } else
     DoeLog(1) && (eLog()<< Verbose(1) << "Atom " << pointer->getName() << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl);
@@ -700 +761 @@
     ElementCount--;
   RemoveBonds(pointer);
-  return remove(pointer, start, end);
+  erase(pointer);
+  return true;
 };
 
@@ -717 +779 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
-  unlink(pointer);
+  erase(pointer);
   return true;
 };
@@ -726 +788 @@
 bool molecule::CleanupMolecule()
 {
-  return (cleanup(first,last) && cleanup(start,end));
+  for (molecule::iterator iter = begin(); !empty(); iter = begin())
+      erase(iter);
 };
 
@@ -733 +796 @@
  * \return pointer to atom or NULL
  */
-atom * molecule::FindAtom(int Nr)  const{
-  atom * walker = find(&Nr, start,end);
-  if (walker != NULL) {
+atom * molecule::FindAtom(int Nr)  const
+{
+  molecule::const_iterator iter = begin();
+  for (; iter != end(); ++iter)
+    if ((*iter)->nr == Nr)
+      break;
+  if (iter != end()) {
     //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
-    return walker;
+    return (*iter);
   } else {
     DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
@@ -867 +934 @@
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     for (int step=0;step<MDSteps;step++) {
-      *output << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
+      *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
       ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     }
@@ -884 +951 @@
   if (output != NULL) {
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
-    *output << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
+    *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
     ActOnAllAtoms( &atom::OutputXYZLine, output );
     return true;
@@ -894 +961 @@
  * \param *out output stream for debugging
  */
-void molecule::CountAtoms()
-{
+int molecule::doCountAtoms()
+{
+  int res = size();
   int i = 0;
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
+  NoNonHydrogen = 0;
+  for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+    (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
+    if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
+      NoNonHydrogen++;
+    stringstream sstr;
+    sstr << (*iter)->type->symbol << (*iter)->nr+1;
+    (*iter)->setName(sstr.str());
+    DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
     i++;
   }
-  if ((AtomCount == 0) || (i != AtomCount)) {
-    DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);
-    AtomCount = i;
-
-    // count NonHydrogen atoms and give each atom a unique name
-    if (AtomCount != 0) {
-      i=0;
-      NoNonHydrogen = 0;
-      Walker = start;
-      while (Walker->next != end) {
-        Walker = Walker->next;
-        Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-        if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
-          NoNonHydrogen++;
-        stringstream sstr;
-        sstr << Walker->type->symbol << Walker->nr+1;
-        Walker->setName(sstr.str());
-        DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->getName() << "." << endl);
-        i++;
-      }
-    } else
-      DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);
-  }
+  return res;
 };
 
@@ -986 +1038 @@
   /// first count both their atoms and elements and update lists thereby ...
   //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
-  CountAtoms();
-  OtherMolecule->CountAtoms();
   CountElements();
   OtherMolecule->CountElements();
@@ -994 +1044 @@
   /// -# AtomCount
   if (result) {
-    if (AtomCount != OtherMolecule->AtomCount) {
-      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
+    if (getAtomCount() != OtherMolecule->getAtomCount()) {
+      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);
       result = false;
-    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+    } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
   }
   /// -# ElementCount
@@ -1034 +1084 @@
   if (result) {
     DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
-    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
-    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
+    Distances = new double[getAtomCount()];
+    OtherDistances = new double[getAtomCount()];
     SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
     SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
+    for(int i=0;i<getAtomCount();i++) {
+      Distances[i] = 0.;
+      OtherDistances[i] = 0.;
+    }
 
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
     DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
-    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
-    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
-    gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
-    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
-    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
+    PermMap = new size_t[getAtomCount()];
+    OtherPermMap = new size_t[getAtomCount()];
+    for(int i=0;i<getAtomCount();i++) {
+      PermMap[i] = 0;
+      OtherPermMap[i] = 0;
+    }
+    gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
+    gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
+    PermutationMap = new int[getAtomCount()];
+    for(int i=0;i<getAtomCount();i++)
+      PermutationMap[i] = 0;
     DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
-    for(int i=AtomCount;i--;)
+    for(int i=getAtomCount();i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
 
@@ -1053 +1113 @@
     DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
     flag = 0;
-    for (int i=0;i<AtomCount;i++) {
+    for (int i=0;i<getAtomCount();i++) {
       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)
@@ -1060 +1120 @@
 
     // free memory
-    Free(&PermMap);
-    Free(&OtherPermMap);
-    Free(&Distances);
-    Free(&OtherDistances);
+    delete[](PermMap);
+    delete[](OtherPermMap);
+    delete[](Distances);
+    delete[](OtherDistances);
     if (flag) { // if not equal
-      Free(&PermutationMap);
+      delete[](PermutationMap);
       result = false;
     }
@@ -1089 +1149 @@
 int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
 {
-  atom *Walker = NULL, *OtherWalker = NULL;
   DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
-  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
-  for (int i=AtomCount;i--;)
+  int *AtomicMap = new int[getAtomCount()];
+  for (int i=getAtomCount();i--;)
     AtomicMap[i] = -1;
   if (OtherMolecule == this) {  // same molecule
-    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
+    for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
       AtomicMap[i] = i;
     DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
   } else {
     DoLog(4) && (Log() << Verbose(4) << "Map is ");
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      if (Walker->father == NULL) {
-        AtomicMap[Walker->nr] = -2;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+      if ((*iter)->father == NULL) {
+        AtomicMap[(*iter)->nr] = -2;
       } else {
-        OtherWalker = OtherMolecule->start;
-        while (OtherWalker->next != OtherMolecule->end) {
-          OtherWalker = OtherWalker->next;
+        for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
       //for (int i=0;i<AtomCount;i++) { // search atom
-        //for (int j=0;j<OtherMolecule->AtomCount;j++) {
-          //Log() << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
-          if (Walker->father == OtherWalker)
-            AtomicMap[Walker->nr] = OtherWalker->nr;
+        //for (int j=0;j<OtherMolecule->getAtomCount();j++) {
+          //Log() << Verbose(4) << "Comparing father " << (*iter)->father << " with the other one " << (*runner)->father << "." << endl;
+          if ((*iter)->father == (*runner))
+            AtomicMap[(*iter)->nr] = (*runner)->nr;
         }
       }
-      DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");
+      DoLog(0) && (Log() << Verbose(0) << AtomicMap[(*iter)->nr] << "\t");
     }
     DoLog(0) && (Log() << Verbose(0) << endl);
@@ -1150 +1205 @@
 void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
 {
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    array[(Walker->*index)] = Walker;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    array[((*iter)->*index)] = (*iter);
   }
 };