/* * 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. */ /* * ConstructBondGraphAction.cpp * * Created on: May 9, 2010 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "Helpers/MemDebug.hpp" #include "Descriptors/AtomIdDescriptor.hpp" #include "Descriptors/MoleculeDescriptor.hpp" #include "atom.hpp" #include "bond.hpp" #include "bondgraph.hpp" #include "config.hpp" #include "linkedcell.hpp" #include "Helpers/Log.hpp" #include "Helpers/Verbose.hpp" #include "molecule.hpp" #include "World.hpp" #include #include typedef std::map< moleculeId_t, std::vector > MolAtomList; using namespace std; #include "Actions/FragmentationAction/ConstructBondGraphAction.hpp" // and construct the stuff #include "ConstructBondGraphAction.def" #include "Action_impl_pre.hpp" /** =========== define the function ====================== */ Action::state_ptr FragmentationConstructBondGraphAction::performCall() { // obtain information getParametersfromValueStorage(); DoLog(1) && (Log() << Verbose(1) << "Constructing bond graph for selected atoms ... " << endl); config *configuration = World::getInstance().getConfig(); BondGraph *BG = configuration->BG; ASSERT(BG != NULL, "FragmentationConstructBondGraphAction: BondGraph is NULL."); double BondDistance = BG->getMaxDistance(); bool IsAngstroem = configuration->GetIsAngstroem(); atom *Walker = NULL; atom *OtherWalker = NULL; int n[NDIM]; double MinDistance, MaxDistance; 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++; for(BondList::iterator BondRunner = (AtomRunner->second)->ListOfBonds.begin(); !(AtomRunner->second)->ListOfBonds.empty(); BondRunner = (AtomRunner->second)->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); LinkedCell::LinkedNodes list; for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection(); AtomRunner != World::getInstance().endAtomSelection(); ++AtomRunner) { list.push_back(AtomRunner->second); } LC = new LinkedCell(list, 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 LinkedCell::LinkedNodes *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 (LinkedCell::LinkedNodes::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 LinkedCell::LinkedNodes *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 (LinkedCell::LinkedNodes::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; BG->CovalentMinMaxDistance(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem); const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition()); const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance); // Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl; 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(Binder); OtherWalker->father->RegisterBond(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 FragmentationConstructBondGraphAction::performUndo(Action::state_ptr _state) { // FragmentationConstructBondGraphState *state = assert_cast(_state.get()); return Action::success; } Action::state_ptr FragmentationConstructBondGraphAction::performRedo(Action::state_ptr _state){ return Action::success; } bool FragmentationConstructBondGraphAction::canUndo() { return false; } bool FragmentationConstructBondGraphAction::shouldUndo() { return false; } const string FragmentationConstructBondGraphAction::getName() { return NAME; } /** =========== end of function ====================== */