Context Navigation

-              r539f32
+              rfa9d1d
 #include "CodePatterns/MemDebug.hpp"
+#include "Descriptors/AtomIdDescriptor.hpp"
+#include "Descriptors/MoleculeDescriptor.hpp"
+#include "CodePatterns/Log.hpp"
+#include "atom.hpp"
+#include "Bond/bond.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Range.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "config.hpp"
+#include "Descriptors/AtomSelectionDescriptor.hpp"
 #include "Graph/BondGraph.hpp"
-#include "linkedcell.hpp"
 #include "molecule.hpp"
-#include "PointCloudAdaptor.hpp"
 #include "World.hpp"
-#include "WorldTime.hpp"
 #include <iostream>
 #include <list>
 #include <string>
-typedef std::map< moleculeId_t, std::vector<atomId_t> > MolAtomList;
 #include "Actions/GraphAction/CreateAdjacencyAction.hpp"
 …
   getParametersfromValueStorage();
-  DoLog(1) && (Log() << Verbose(1) << "Constructing bond graph for selected atoms ... " << endl);
   BondGraph *BG = World::getInstance().getBondGraph();
   ASSERT(BG != NULL, "GraphCreateAdjacencyAction: BondGraph is NULL.");
-  double BondDistance = BG->getMaxPossibleBondDistance(AtomSetMixin<std::vector<atom *> >(World::getInstance().getSelectedAtoms()));
+  atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  int n[NDIM];
+  LinkedCell *LC = NULL;
+  Box &domain = World::getInstance().getDomain();
+  World::AtomComposite Set = World::getInstance().getAllAtoms(AtomsBySelection());
+  BG->CreateAdjacency(Set);
+  // remove every bond from the selected atoms' list
+  int AtomCount = 0;
+  for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
+      AtomRunner != World::getInstance().endAtomSelection();
+      ++AtomRunner) {
+    AtomCount++;
+    (AtomRunner->second)->removeAllBonds();
+//    BondList& ListOfBonds = (AtomRunner->second)->getListOfBonds();
+//    for(BondList::iterator BondRunner = ListOfBonds.begin();
+//        !ListOfBonds.empty();
+//        BondRunner = ListOfBonds.begin())
+//      if ((*BondRunner)->leftatom == AtomRunner->second)
+//        delete((*BondRunner));
+  }
+  int BondCount = 0;
+  // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
+  DoLog(1) && (Log() << Verbose(1) << "AtomCount " << AtomCount << " and bonddistance is " << BondDistance << "." << endl);
+  if ((AtomCount > 1) && (BondDistance > 1.)) {
+    DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
+    TesselPointSTLList list;
+    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
+        AtomRunner != World::getInstance().endAtomSelection();
+        ++AtomRunner) {
+      list.push_back(AtomRunner->second);
+    }
+    PointCloudAdaptor< TesselPointSTLList > cloud(&list, "AtomSelection");
+    LC = new LinkedCell(cloud, BondDistance);
+    // create a list to map Tesselpoint::Nr to atom *
+    DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl);
+    // set numbers for atoms that can later be used
+    std::map<TesselPoint *, int> AtomIds;
+    int i=0;
+    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
+        AtomRunner != World::getInstance().endAtomSelection();
+        ++AtomRunner) {
+      AtomIds.insert(pair<TesselPoint *, int> (AtomRunner->second, i++) );
+    }
+    // 3a. go through every cell
+    DoLog(2) && (Log() << Verbose(2) << "Celling ... " << endl);
+    for (LC->n[0] = 0; LC->n[0] < LC->N[0]; LC->n[0]++)
+      for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
+        for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
+          const TesselPointSTLList *List = LC->GetCurrentCell();
+//          Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+          if (List != NULL) {
+            for (TesselPointSTLList::const_iterator Runner = List->begin();
+                Runner != List->end();
+                Runner++) {
+              Walker = dynamic_cast<atom*>(*Runner);
+              ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
+              //Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
+              // 3c. check for possible bond between each atom in this and every one in the 27 cells
+              for (n[0] = -1; n[0] <= 1; n[0]++)
+                for (n[1] = -1; n[1] <= 1; n[1]++)
+                  for (n[2] = -1; n[2] <= 1; n[2]++) {
+                    const TesselPointSTLList *OtherList = LC->GetRelativeToCurrentCell(n);
+//                    Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+                    if (OtherList != NULL) {
+                      for (TesselPointSTLList::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
+                        if (AtomIds.find(*OtherRunner)->second > AtomIds.find(Walker)->second) {
+                          OtherWalker = dynamic_cast<atom*>(*OtherRunner);
+                          ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
+                          //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
+                          const range<double> MinMaxDistanceSquared(
+                              BG->getMinMaxDistanceSquared(Walker, OtherWalker));
+                          const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
+                          const bool status = MinMaxDistanceSquared.isInRange(distance);
+//                          LOG3, "INFO: MinMaxDistanceSquared interval is " << MinMaxDistanceSquared << ".");
+                          if (AtomIds[OtherWalker->father] > AtomIds[Walker->father]) {
+                            if (status) { // create bond if distance is smaller
+//                              Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
+                              //const bond * Binder =
+                                  Walker->father->addBond(WorldTime::getTime(), OtherWalker->father);
+                              BondCount++;
+                            } else {
+//                              Log() << Verbose(1) << "Not Adding: distance too great." << endl;
+                            }
+                          } else {
+//                            Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
+                          }
+                        }
+                      }
+                    }
+                  }
+            }
+          }
+        }
+    delete (LC);
+    DoLog(1) && (Log() << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl);
+    // correct bond degree by comparing valence and bond degree
+    DoLog(2) && (Log() << Verbose(2) << "Correcting bond degree ... " << endl);
+    //CorrectBondDegree();
+  } else
+    DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl);
+  DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl);
+  size_t BondCount = 0;
+  std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
+  for (std::vector<molecule *>::const_iterator iter = molecules.begin();
+      iter != molecules.end(); ++iter)
+    BondCount += (*iter)->getBondCount();
+  LOG(0, "STATUS: Recognized " << BondCount << " bonds.");
   return Action::success;

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset fa9d1d for src/Actions/GraphAction/CreateAdjacencyAction.cpp

Legend:

src/Actions/GraphAction/CreateAdjacencyAction.cpp

Download in other formats: