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

File:

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

src/molecule.cpp

@@ -4 +4 @@
  *
  */
-
-#include "Helpers/MemDebug.hpp"
 
 #include <cstring>
@@ -37 +35 @@
  * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
  */
-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())
-{
+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);
 
   // other stuff
@@ -61 +67 @@
 {
   CleanupMolecule();
+  delete(first);
+  delete(last);
+  end->getWorld()->destroyAtom(end);
+  start->getWorld()->destroyAtom(start);
 };
 
@@ -73 +83 @@
 }
 
-int molecule::getAtomCount() const{
-  return *AtomCount;
-}
-
 void molecule::setName(const std::string _name){
   OBSERVE;
@@ -98 +104 @@
   stringstream sstr;
   periodentafel *periode = World::getInstance().getPeriode();
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    counts[(*iter)->type->getNumber()]++;
+  for(atom *Walker = start; Walker != end; Walker = Walker->next) {
+    counts[Walker->type->getNumber()]++;
   }
   std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
@@ -109 +115 @@
 }
 
-/************************** 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.
@@ -184 +123 @@
 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)
@@ -199 +141 @@
       }
     }
-    insert(pointer);
-    pointer->setMolecule(this);
-  }
-  return true;
+    retval = add(pointer, end);
+  }
+  return retval;
 };
 
@@ -216 +157 @@
   if (pointer != NULL) {
     atom *walker = pointer->clone();
-    walker->setName(pointer->getName());
+    stringstream sstr;
+    sstr << pointer->getName();
+    walker->setName(sstr.str());
     walker->nr = last_atom++;  // increase number within molecule
-    insert(walker);
+    add(walker, end);
     if ((pointer->type != NULL) && (pointer->type->Z != 1))
       NoNonHydrogen++;
+    AtomCount++;
     retval=walker;
   }
@@ -298 +242 @@
     Orthovector1.MatrixMultiplication(matrix);
     InBondvector -= Orthovector1; // subtract just the additional translation
-    delete[](matrix);
+    Free(&matrix);
     bondlength = InBondvector.Norm();
 //    Log() << Verbose(4) << "Corrected InBondvector is now: ";
@@ -307 +251 @@
   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) {
@@ -529 +472 @@
       break;
   }
-  delete[](matrix);
+  Free(&matrix);
 
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
@@ -612 +555 @@
 
   // copy all bonds
-  bond *Binder = NULL;
+  bond *Binder = first;
   bond *NewBond = NULL;
-  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;
-      }
+  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;
+  }
   // correct fathers
   ActOnAllAtoms( &atom::CorrectFather );
 
   // copy values
+  copy->CountAtoms();
   copy->CountElements();
-  if (hasBondStructure()) {  // if adjaceny list is present
+  if (first->next != last) {  // if adjaceny list is present
     copy->BondDistance = BondDistance;
   }
@@ -666 +608 @@
 bond * molecule::AddBond(atom *atom1, atom *atom2, int degree)
 {
-  OBSERVE;
   bond *Binder = NULL;
-
-  // 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++;
-
+  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);
+  }
   return Binder;
 };
@@ -692 +630 @@
 {
   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
-  delete(pointer);
+  pointer->leftatom->RegisterBond(pointer);
+  pointer->rightatom->RegisterBond(pointer);
+  removewithoutcheck(pointer);
   return true;
 };
@@ -752 +692 @@
 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);
@@ -761 +700 @@
     ElementCount--;
   RemoveBonds(pointer);
-  erase(pointer);
-  return true;
+  return remove(pointer, start, end);
 };
 
@@ -779 +717 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
-  erase(pointer);
+  unlink(pointer);
   return true;
 };
@@ -788 +726 @@
 bool molecule::CleanupMolecule()
 {
-  for (molecule::iterator iter = begin(); !empty(); iter = begin())
-      erase(iter);
+  return (cleanup(first,last) && cleanup(start,end));
 };
 
@@ -796 +733 @@
  * \return pointer to atom or NULL
  */
-atom * molecule::FindAtom(int Nr)  const
-{
-  molecule::const_iterator iter = begin();
-  for (; iter != end(); ++iter)
-    if ((*iter)->nr == Nr)
-      break;
-  if (iter != end()) {
+atom * molecule::FindAtom(int Nr)  const{
+  atom * walker = find(&Nr, start,end);
+  if (walker != NULL) {
     //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
-    return (*iter);
+    return walker;
   } else {
     DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
@@ -934 +867 @@
     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 << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
+      *output << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
       ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     }
@@ -951 +884 @@
   if (output != NULL) {
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
-    *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
+    *output << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
     ActOnAllAtoms( &atom::OutputXYZLine, output );
     return true;
@@ -961 +894 @@
  * \param *out output stream for debugging
  */
-int molecule::doCountAtoms()
-{
-  int res = size();
+void molecule::CountAtoms()
+{
   int i = 0;
-  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);
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
     i++;
   }
-  return res;
+  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);
+  }
 };
 
@@ -1038 +986 @@
   /// 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();
@@ -1044 +994 @@
   /// -# AtomCount
   if (result) {
-    if (getAtomCount() != OtherMolecule->getAtomCount()) {
-      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);
+    if (AtomCount != OtherMolecule->AtomCount) {
+      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
       result = false;
-    } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
+    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
   }
   /// -# ElementCount
@@ -1084 +1034 @@
   if (result) {
     DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
-    Distances = new double[getAtomCount()];
-    OtherDistances = new double[getAtomCount()];
+    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
+    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
     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 = 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;
+    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");
     DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
-    for(int i=getAtomCount();i--;)
+    for(int i=AtomCount;i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
 
@@ -1113 +1053 @@
     DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
     flag = 0;
-    for (int i=0;i<getAtomCount();i++) {
+    for (int i=0;i<AtomCount;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)
@@ -1120 +1060 @@
 
     // free memory
-    delete[](PermMap);
-    delete[](OtherPermMap);
-    delete[](Distances);
-    delete[](OtherDistances);
+    Free(&PermMap);
+    Free(&OtherPermMap);
+    Free(&Distances);
+    Free(&OtherDistances);
     if (flag) { // if not equal
-      delete[](PermutationMap);
+      Free(&PermutationMap);
       result = false;
     }
@@ -1149 +1089 @@
 int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
 {
+  atom *Walker = NULL, *OtherWalker = NULL;
   DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
-  int *AtomicMap = new int[getAtomCount()];
-  for (int i=getAtomCount();i--;)
+  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
+  for (int i=AtomCount;i--;)
     AtomicMap[i] = -1;
   if (OtherMolecule == this) {  // same molecule
-    for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
+    for (int i=AtomCount;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 ");
-    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-      if ((*iter)->father == NULL) {
-        AtomicMap[(*iter)->nr] = -2;
+    Walker = start;
+    while (Walker->next != end) {
+      Walker = Walker->next;
+      if (Walker->father == NULL) {
+        AtomicMap[Walker->nr] = -2;
       } else {
-        for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
+        OtherWalker = OtherMolecule->start;
+        while (OtherWalker->next != OtherMolecule->end) {
+          OtherWalker = OtherWalker->next;
       //for (int i=0;i<AtomCount;i++) { // search atom
-        //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;
+        //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;
         }
       }
-      DoLog(0) && (Log() << Verbose(0) << AtomicMap[(*iter)->nr] << "\t");
+      DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");
     }
     DoLog(0) && (Log() << Verbose(0) << endl);
@@ -1205 +1150 @@
 void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
 {
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    array[((*iter)->*index)] = (*iter);
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    array[(Walker->*index)] = Walker;
   }
 };