/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * CreateAdjacencyAction.cpp * * Created on: May 9, 2010 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "Descriptors/AtomIdDescriptor.hpp" #include "Descriptors/MoleculeDescriptor.hpp" #include "atom.hpp" #include "Bond/bond.hpp" #include "CodePatterns/Log.hpp" #include "CodePatterns/Range.hpp" #include "CodePatterns/Verbose.hpp" #include "config.hpp" #include "Graph/BondGraph.hpp" #include "linkedcell.hpp" #include "molecule.hpp" #include "PointCloudAdaptor.hpp" #include "World.hpp" #include "WorldTime.hpp" #include #include #include typedef std::map< moleculeId_t, std::vector > MolAtomList; using namespace std; #include "Actions/FragmentationAction/CreateAdjacencyAction.hpp" // and construct the stuff #include "CreateAdjacencyAction.def" #include "Action_impl_pre.hpp" /** =========== define the function ====================== */ Action::state_ptr FragmentationCreateAdjacencyAction::performCall() { // obtain information getParametersfromValueStorage(); DoLog(1) && (Log() << Verbose(1) << "Constructing bond graph for selected atoms ... " << endl); BondGraph *BG = World::getInstance().getBondGraph(); ASSERT(BG != NULL, "FragmentationCreateAdjacencyAction: BondGraph is NULL."); double BondDistance = BG->getMaxPossibleBondDistance(AtomSetMixin >(World::getInstance().getSelectedAtoms())); atom *Walker = NULL; atom *OtherWalker = NULL; int n[NDIM]; LinkedCell *LC = NULL; Box &domain = World::getInstance().getDomain(); // remove every bond from the selected atoms' list int AtomCount = 0; for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection(); AtomRunner != World::getInstance().endAtomSelection(); ++AtomRunner) { AtomCount++; 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 AtomIds; int i=0; for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection(); AtomRunner != World::getInstance().endAtomSelection(); ++AtomRunner) { AtomIds.insert(pair (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(*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(*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 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; bond * const Binder = new bond(Walker->father, OtherWalker->father, 1); BondCount++; Walker->father->RegisterBond(WorldTime::getTime(),Binder); OtherWalker->father->RegisterBond(WorldTime::getTime(),Binder); } 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); return Action::success; } Action::state_ptr FragmentationCreateAdjacencyAction::performUndo(Action::state_ptr _state) { // FragmentationCreateAdjacencyState *state = assert_cast(_state.get()); return Action::success; } Action::state_ptr FragmentationCreateAdjacencyAction::performRedo(Action::state_ptr _state){ return Action::success; } bool FragmentationCreateAdjacencyAction::canUndo() { return false; } bool FragmentationCreateAdjacencyAction::shouldUndo() { return false; } /** =========== end of function ====================== */