Changes in src/molecule.cpp [952f38:822f01]

File:

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

src/molecule.cpp

@@ -24 +24 @@
 #include "element.hpp"
 #include "graph.hpp"
-#include "Helpers/helpers.hpp"
+#include "helpers.hpp"
 #include "leastsquaremin.hpp"
 #include "linkedcell.hpp"
 #include "lists.hpp"
-#include "Helpers/Log.hpp"
+#include "log.hpp"
 #include "molecule.hpp"
 
@@ -34 +34 @@
 #include "stackclass.hpp"
 #include "tesselation.hpp"
-#include "LinearAlgebra/Vector.hpp"
-#include "LinearAlgebra/Matrix.hpp"
+#include "vector.hpp"
+#include "Matrix.hpp"
 #include "World.hpp"
 #include "Box.hpp"
-#include "LinearAlgebra/Plane.hpp"
+#include "Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
 
@@ -151 +151 @@
 molecule::const_iterator molecule::erase( const_iterator loc )
 {
+  OBSERVE;
   molecule::const_iterator iter = loc;
   iter--;
@@ -156 +157 @@
   atomIds.erase( atom->getId() );
   atoms.remove( atom );
+  formula-=atom->type;
   atom->removeFromMolecule();
   return iter;
@@ -162 +164 @@
 molecule::const_iterator molecule::erase( atom * key )
 {
+  OBSERVE;
   molecule::const_iterator iter = find(key);
   if (iter != end()){
     atomIds.erase( key->getId() );
     atoms.remove( key );
+    formula-=key->type;
     key->removeFromMolecule();
   }
@@ -183 +187 @@
 pair<molecule::iterator,bool> molecule::insert ( atom * const key )
 {
+  OBSERVE;
   pair<atomIdSet::iterator,bool> res = atomIds.insert(key->getId());
   if (res.second) { // push atom if went well
     atoms.push_back(key);
+    formula+=key->type;
     return pair<iterator,bool>(molecule::iterator(--end()),res.second);
   } else {
@@ -212 +218 @@
       if(pointer->getName() == "Unknown"){
         stringstream sstr;
-        sstr << pointer->type->symbol << pointer->nr+1;
+        sstr << pointer->type->getSymbol() << pointer->nr+1;
         pointer->setName(sstr.str());
       }
@@ -233 +239 @@
   if (pointer != NULL) {
     atom *walker = pointer->clone();
+    formula += walker->type;
     walker->setName(pointer->getName());
     walker->nr = last_atom++;  // increase number within molecule
@@ -619 +626 @@
 {
   molecule *copy = World::getInstance().createMolecule();
-  atom *LeftAtom = NULL, *RightAtom = NULL;
 
   // copy all atoms
-  ActOnCopyWithEachAtom ( &molecule::AddCopyAtom, copy );
+  for_each(atoms.begin(),atoms.end(),bind1st(mem_fun(&molecule::AddCopyAtom),copy));
 
   // copy all bonds
-  bond *Binder = NULL;
-  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);
+        bond *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);
+        atomSet::iterator leftiter=find_if(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->leftatom));
+        atomSet::iterator rightiter=find_if(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->rightatom));
+        ASSERT(leftiter!=atoms.end(),"No original left atom for bondcopy found");
+        ASSERT(leftiter!=atoms.end(),"No original right atom for bondcopy found");
+        atom *LeftAtom = *leftiter;
+        atom *RightAtom = *rightiter;
+
+        bond *NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
         NewBond->Cyclic = Binder->Cyclic;
         if (Binder->Cyclic)
@@ -643 +651 @@
       }
   // correct fathers
-  ActOnAllAtoms( &atom::CorrectFather );
+  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
 
   // copy values
@@ -852 +860 @@
  * \param *out output stream
  */
-bool molecule::Output(ofstream * const output)
-{
-  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
-
-  for (int i=0;i<MAX_ELEMENTS;++i) {
-    AtomNo[i] = 0;
-    ElementNo[i] = 0;
-  }
+bool molecule::Output(ostream * const output)
+{
   if (output == NULL) {
     return false;
   } else {
+    int AtomNo[MAX_ELEMENTS];
+    memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
+    enumeration<const element*> elementLookup = formula.enumerateElements();
+    for(map<const element*,unsigned int>::iterator iter=elementLookup.there.begin();
+        iter!=elementLookup.there.end();++iter){
+      cout << "Enumerated element " << *iter->first << " with number " << iter->second << 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;
-    for (int i=0;i<MAX_ELEMENTS;++i) {
-      if (ElementNo[i] == 1)
-        ElementNo[i] = current++;
-    }
-    ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
+    for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
     return true;
   }
@@ -913 +916 @@
 {
   DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
-  ActOnAllAtoms (&atom::OutputBondOfAtom );
+  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
   DoLog(0) && (Log() << Verbose(0) << endl);
 };
@@ -936 +939 @@
     for (int step=0;step<MDSteps;step++) {
       *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
-      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
+      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
     }
     return true;
@@ -953 +956 @@
     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);
-    ActOnAllAtoms( &atom::OutputXYZLine, output );
+    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
     return true;
   } else
@@ -972 +975 @@
       NoNonHydrogen++;
     stringstream sstr;
-    sstr << (*iter)->type->symbol << (*iter)->nr+1;
+    sstr << (*iter)->type->getSymbol() << (*iter)->nr+1;
     (*iter)->setName(sstr.str());
     DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
@@ -978 +981 @@
   }
   return res;
-};
-
-/** Determines whether two molecules actually contain the same atoms and coordination.
- * \param *out output stream for debugging
- * \param *OtherMolecule the molecule to compare this one to
- * \param threshold upper limit of difference when comparing the coordination.
- * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
- */
-int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)
-{
-  int flag;
-  double *Distances = NULL, *OtherDistances = NULL;
-  Vector CenterOfGravity, OtherCenterOfGravity;
-  size_t *PermMap = NULL, *OtherPermMap = NULL;
-  int *PermutationMap = NULL;
-  bool result = true; // status of comparison
-
-  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;
-
-  /// ... and compare:
-  /// -# AtomCount
-  if (result) {
-    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: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
-  }
-  /// -# Formula
-  if (result) {
-    if (formula != OtherMolecule->formula) {
-      DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);
-      result = false;
-    } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;
-  }
-  /// then determine and compare center of gravity for each molecule ...
-  if (result) {
-    DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
-    DeterminePeriodicCenter(CenterOfGravity);
-    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
-    DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl);
-    DoLog(5) && (Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl);
-    if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {
-      DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
-      result = false;
-    }
-  }
-
-  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
-  if (result) {
-    DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
-    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 = 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=getAtomCount();i--;)
-      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
-
-    /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
-    DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
-    flag = 0;
-    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)
-        flag = 1;
-    }
-
-    // free memory
-    delete[](PermMap);
-    delete[](OtherPermMap);
-    delete[](Distances);
-    delete[](OtherDistances);
-    if (flag) { // if not equal
-      delete[](PermutationMap);
-      result = false;
-    }
-  }
-  /// return pointer to map if all distances were below \a threshold
-  DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
-  if (result) {
-    DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
-    return PermutationMap;
-  } else {
-    DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
-    return NULL;
-  }
 };