Changes in src/molecule.cpp [f8e486:7baf4a]

File:

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

src/molecule.cpp

@@ -5 +5 @@
  */
 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
 #include "Helpers/MemDebug.hpp"
 
 #include <cstring>
 #include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_heapsort.h>
 
 #include "World.hpp"
@@ -22 +30 @@
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
+
 #include "periodentafel.hpp"
 #include "stackclass.hpp"
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 #include "Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
@@ -39 +49 @@
 molecule::molecule(const periodentafel * const teil) :
   Observable("molecule"),
-  elemente(teil),  MDSteps(0),  BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0),
+  elemente(teil),  MDSteps(0),  BondCount(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
-  for(int i=MAX_ELEMENTS;i--;)
-    ElementsInMolecule[i] = 0;
-  strcpy(name,World::getInstance().getDefaultName());
+  AtomCount(this,boost::bind(&molecule::doCountAtoms,this),"AtomCount"), last_atom(0),  InternalPointer(atoms.begin())
+{
+
+  strcpy(name,World::getInstance().getDefaultName().c_str());
 };
 
@@ -79 +85 @@
 void molecule::setName(const std::string _name){
   OBSERVE;
+  cout << "Set name of molecule " << getId() << " to " << _name << endl;
   strncpy(name,_name.c_str(),MAXSTRINGSIZE);
 }
@@ -90 +97 @@
 }
 
-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 (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    counts[(*iter)->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();
+const Formula &molecule::getFormula(){
+  return formula;
+}
+
+unsigned int molecule::getElementCount(){
+  return formula.getElementCount();
+}
+
+bool molecule::hasElement(const element *element) const{
+  return formula.hasElement(element);
+}
+
+bool molecule::hasElement(atomicNumber_t Z) const{
+  return formula.hasElement(Z);
+}
+
+bool molecule::hasElement(const string &shorthand) const{
+  return formula.hasElement(shorthand);
 }
 
@@ -143 +151 @@
 molecule::const_iterator molecule::erase( const_iterator loc )
 {
+  OBSERVE;
   molecule::const_iterator iter = loc;
   iter--;
   atom* atom = *loc;
-  atoms.erase( loc );
+  atomIds.erase( atom->getId() );
+  atoms.remove( atom );
+  formula-=atom->type;
   atom->removeFromMolecule();
   return iter;
@@ -153 +164 @@
 molecule::const_iterator molecule::erase( atom * key )
 {
-  cout << "trying to erase atom" << endl;
+  OBSERVE;
   molecule::const_iterator iter = find(key);
   if (iter != end()){
-    atoms.erase( iter++ );
+    atomIds.erase( key->getId() );
+    atoms.remove( key );
+    formula-=key->type;
     key->removeFromMolecule();
   }
@@ -164 +177 @@
 molecule::const_iterator molecule::find ( atom * key ) const
 {
-  return atoms.find( key );
+  molecule::const_iterator iter;
+  for (molecule::const_iterator Runner = begin(); Runner != end(); ++Runner) {
+    if (*Runner == key)
+      return molecule::const_iterator(Runner);
+  }
+  return molecule::const_iterator(atoms.end());
 }
 
 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);
+  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 {
+    return pair<iterator,bool>(molecule::iterator(end()),res.second);
+  }
 }
 
 bool molecule::containsAtom(atom* key){
-  return atoms.count(key);
+  return (find(key) != end());
 }
 
@@ -188 +213 @@
     pointer->sort = &pointer->nr;
     if (pointer->type != NULL) {
-      if (ElementsInMolecule[pointer->type->Z] == 0)
-        ElementCount++;
-      ElementsInMolecule[pointer->type->Z]++; // increase number of elements
+      formula += pointer->type;
       if (pointer->type->Z != 1)
         NoNonHydrogen++;
@@ -216 +239 @@
   if (pointer != NULL) {
     atom *walker = pointer->clone();
+    formula += walker->type;
     walker->setName(pointer->getName());
     walker->nr = last_atom++;  // increase number within molecule
@@ -272 +296 @@
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
-  double *matrix = NULL;
+  const Matrix &matrix =  World::getInstance().getDomain().getM();
   bond *Binder = NULL;
-  double * const cell_size = World::getInstance().getDomain();
 
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -295 +318 @@
       } // (signs are correct, was tested!)
     }
-    matrix = ReturnFullMatrixforSymmetric(cell_size);
-    Orthovector1.MatrixMultiplication(matrix);
+    Orthovector1 *= matrix;
     InBondvector -= Orthovector1; // subtract just the additional translation
-    delete[](matrix);
     bondlength = InBondvector.Norm();
 //    Log() << Verbose(4) << "Corrected InBondvector is now: ";
@@ -529 +550 @@
       break;
   }
-  delete[](matrix);
 
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
@@ -606 +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)
@@ -630 +651 @@
       }
   // correct fathers
-  ActOnAllAtoms( &atom::CorrectFather );
+  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
 
   // copy values
-  copy->CountElements();
   if (hasBondStructure()) {  // if adjaceny list is present
     copy->BondDistance = BondDistance;
@@ -648 +668 @@
  * @param three vectors forming the matrix that defines the shape of the parallelpiped
  */
-molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
+molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
   molecule *copy = World::getInstance().createMolecule();
 
-  ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
+  BOOST_FOREACH(atom *iter,atoms){
+    if(iter->IsInShape(region)){
+      copy->AddCopyAtom(iter);
+    }
+  }
 
   //TODO: copy->BuildInducedSubgraph(this);
@@ -728 +752 @@
   else
     length = strlen(molname) - strlen(endname);
+  cout << "Set name of molecule " << getId() << " to " << molname << endl;
   strncpy(name, molname, length);
   name[length]='\0';
@@ -737 +762 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
-  double * const cell_size = World::getInstance().getDomain();
-  cell_size[0] = dim->at(0);
-  cell_size[1] = 0.;
-  cell_size[2] = dim->at(1);
-  cell_size[3] = 0.;
-  cell_size[4] = 0.;
-  cell_size[5] = dim->at(2);
+  Matrix domain;
+  for(int i =0; i<NDIM;++i)
+    domain.at(i,i) = dim->at(i);
+  World::getInstance().setDomain(domain);
 };
 
@@ -754 +776 @@
   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
-  } 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);
-  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
-    ElementCount--;
+  formula-=pointer->type;
   RemoveBonds(pointer);
   erase(pointer);
@@ -773 +790 @@
   if (pointer == NULL)
     return false;
-  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
-    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
-  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);
-  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
-    ElementCount--;
+  formula-=pointer->type;
   erase(pointer);
   return true;
@@ -790 +802 @@
   for (molecule::iterator iter = begin(); !empty(); iter = begin())
       erase(iter);
+  return empty();
 };
 
@@ -835 +848 @@
 bool molecule::CheckBounds(const Vector *x) const
 {
-  double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   bool result = true;
-  int j =-1;
   for (int i=0;i<NDIM;i++) {
-    j += i+1;
-    result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
+    result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
   }
   //return result;
@@ -849 +860 @@
  * \param *out output stream
  */
-bool molecule::Output(ofstream * const output)
-{
-  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
-  CountElements();
-
-  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, 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;
   }
@@ -879 +884 @@
 {
   int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
-  CountElements();
 
   if (output == NULL) {
@@ -912 +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);
 };
@@ -921 +925 @@
 bool molecule::Checkout(ofstream * const output)  const
 {
-  return elemente->Checkout(output, ElementsInMolecule);
+  return formula.checkOut(output);
 };
 
@@ -935 +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;
@@ -952 +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
@@ -979 +983 @@
 };
 
-/** Brings molecule::ElementCount and molecule::ElementsInMolecule up-to-date.
- */
-void molecule::CountElements()
-{
-  for(int i=MAX_ELEMENTS;i--;)
-    ElementsInMolecule[i] = 0;
-  ElementCount = 0;
-
-  SetIndexedArrayForEachAtomTo ( ElementsInMolecule, &element::Z, &Increment, 1);
-
-  for(int i=MAX_ELEMENTS;i--;)
-    ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
-};
-
-
-/** Counts necessary number of valence electrons and returns number and SpinType.
- * \param configuration containing everything
- */
-void molecule::CalculateOrbitals(class config &configuration)
-{
-  configuration.MaxPsiDouble = configuration.PsiMaxNoDown = configuration.PsiMaxNoUp = configuration.PsiType = 0;
-  for(int i=MAX_ELEMENTS;i--;) {
-    if (ElementsInMolecule[i] != 0) {
-      //Log() << Verbose(0) << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
-      configuration.MaxPsiDouble += ElementsInMolecule[i]*((int)elemente->FindElement(i)->Valence);
-    }
-  }
-  configuration.PsiMaxNoDown = configuration.MaxPsiDouble/2 + (configuration.MaxPsiDouble % 2);
-  configuration.PsiMaxNoUp = configuration.MaxPsiDouble/2;
-  configuration.MaxPsiDouble /= 2;
-  configuration.PsiType = (configuration.PsiMaxNoDown == configuration.PsiMaxNoUp) ? 0 : 1;
-  if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2)) {
-    configuration.ProcPEGamma /= 2;
-    configuration.ProcPEPsi *= 2;
-  } else {
-    configuration.ProcPEGamma *= configuration.ProcPEPsi;
-    configuration.ProcPEPsi = 1;
-  }
-  configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
-};
-
 /** Determines whether two molecules actually contain the same atoms and coordination.
  * \param *out output stream for debugging
@@ -1038 +1001 @@
   /// first count both their atoms and elements and update lists thereby ...
   //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
-  CountElements();
-  OtherMolecule->CountElements();
 
   /// ... and compare:
@@ -1049 +1010 @@
     } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
   }
-  /// -# ElementCount
+  /// -# Formula
   if (result) {
-    if (ElementCount != OtherMolecule->ElementCount) {
-      DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl);
+    if (formula != OtherMolecule->formula) {
+      DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);
       result = false;
-    } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
-  }
-  /// -# ElementsInMolecule
-  if (result) {
-    for (flag=MAX_ELEMENTS;flag--;) {
-      //Log() << Verbose(5) << "Element " <<  flag << ": " << ElementsInMolecule[flag] << " <-> " << OtherMolecule->ElementsInMolecule[flag] << "." << endl;
-      if (ElementsInMolecule[flag] != OtherMolecule->ElementsInMolecule[flag])
-        break;
-    }
-    if (flag < MAX_ELEMENTS) {
-      DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl);
-      result = false;
-    } else Log() << Verbose(4) << "ElementsInMolecule match." << endl;
+    } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;
   }
   /// then determine and compare center of gravity for each molecule ...