/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010-2012 University of Bonn. All rights reserved. * Copyright (C) 2013 Frederik Heber. All rights reserved. * * * This file is part of MoleCuilder. * * MoleCuilder is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * MoleCuilder is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MoleCuilder. If not, see . */ /* * molecule_graph.cpp * * Created on: Oct 5, 2009 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif //#include "CodePatterns/MemDebug.hpp" #include #include "Atom/atom.hpp" #include "Bond/bond.hpp" #include "Box.hpp" #include "CodePatterns/Assert.hpp" #include "CodePatterns/Info.hpp" #include "CodePatterns/Log.hpp" #include "CodePatterns/Verbose.hpp" #include "config.hpp" #include "Graph/DepthFirstSearchAnalysis.hpp" #include "Element/element.hpp" #include "Graph/BondGraph.hpp" #include "Graph/ListOfLocalAtoms.hpp" #include "Helpers/defs.hpp" #include "Helpers/helpers.hpp" #include "LinearAlgebra/RealSpaceMatrix.hpp" #include "LinkedCell/linkedcell.hpp" #include "LinkedCell/PointCloudAdaptor.hpp" #include "molecule.hpp" #include "World.hpp" #include "WorldTime.hpp" /** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule. * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL. * \param *reference reference molecule with the bond structure to be copied * \param ListOfLocalAtoms Lookup table for this subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled * \param FreeList true - ListOfLocalAtoms is free'd before return, false - it is not * \return true - success, false - failure */ bool molecule::FillBondStructureFromReference(const molecule * const reference, ListOfLocalAtoms_t &ListOfLocalAtoms, bool FreeList) { bool status = true; LOG(1, "Begin of FillBondStructureFromReference."); // fill ListOfLocalAtoms if NULL was given if (!FillListOfLocalAtoms(ListOfLocalAtoms, reference->getAtomCount())) { LOG(1, "Filling of ListOfLocalAtoms failed."); return false; } if (status) { LOG(1, "Creating adjacency list for molecule " << getName() << "."); // remove every bond from the list for_each(begin(), end(), boost::bind(&BondedParticle::ClearBondsAtStep,_1,WorldTime::getTime())); for(molecule::iterator iter = begin(); iter != end(); ++iter) { const atom * const Father = (*iter)->GetTrueFather(); //const int AtomNo = Father->getNr(); // global id of the current walker const BondList& ListOfBonds = Father->getListOfBonds(); for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) { atom * const OtherAtom = (*Runner)->GetOtherAtom((*iter)->GetTrueFather()); const ListOfLocalAtoms_t::const_iterator localiter = ListOfLocalAtoms.find(OtherAtom->getNr()); ASSERT( localiter != ListOfLocalAtoms.end(), "molecule::FillBondStructureFromReference() - could not find id" +toString(OtherAtom->getNr())+" in ListOfLocalAtoms."); atom * const OtherWalker = localiter->second; // local copy of current bond partner of walker if (OtherWalker != NULL) { if (OtherWalker->getNr() > (*iter)->getNr()) AddBond((*iter), OtherWalker, (*Runner)->getDegree()); } else { LOG(1, "OtherWalker = ListOfLocalAtoms[" << OtherAtom->getNr() << "] is NULL!"); status = false; } } } } if ((FreeList) && (!ListOfLocalAtoms.empty())) { // free the index lookup list ListOfLocalAtoms.clear(); } LOG(1, "End of FillBondStructureFromReference."); return status; }; /** Checks for presence of bonds within atom list. * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...) * \return true - bonds present, false - no bonds */ bool molecule::hasBondStructure() const { for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) { //LOG(0, "molecule::hasBondStructure() - checking bond list of atom " << (*AtomRunner)->getId() << "."); const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds(); if (!ListOfBonds.empty()) return true; } return false; } /** Prints a list of all bonds to \a *out. */ void molecule::OutputBondsList() const { if (DoLog(1)) { std::stringstream output; output << std::endl << "From contents of bond chain list:"; for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) { const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds(); for(BondList::const_iterator BondRunner = ListOfBonds.begin(); BondRunner != ListOfBonds.end(); ++BondRunner) if ((*BondRunner)->leftatom == *AtomRunner) { output << *(*BondRunner) << "\t"; } } LOG(1, output.str()); } } /** Storing the bond structure of a molecule to file. * Simply stores Atom::Nr and then the Atom::Nr of all bond partners, one per line. * \param &filename name of file * \param path path to file, defaults to empty * \return true - file written successfully, false - writing failed */ bool molecule::StoreBondsToFile(std::string filename, std::string path) { ofstream BondFile; string line; bool status = true; if (path != "") line = path + "/" + filename; else line = filename; BondFile.open(line.c_str(), ios::out); LOG(1, "Saving adjacency list ... "); if (BondFile.good()) { BondFile << "m\tn" << endl; for_each(begin(),end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile)); BondFile.close(); LOG(1, "\t... done."); } else { LOG(1, "\t... failed to open file " << line << "."); status = false; } return status; } ; /** Adds a bond as a copy to a given one * \param *left leftatom of new bond * \param *right rightatom of new bond * \param *CopyBond rest of fields in bond are copied from this * \return pointer to new bond */ bond::ptr molecule::CopyBond(atom *left, atom *right, bond::ptr CopyBond) { bond::ptr Binder = AddBond(left, right, CopyBond->getDegree()); Binder->Cyclic = CopyBond->Cyclic; Binder->Type = CopyBond->Type; return Binder; } ; /** Fills a lookup list of father's Atom::nr -> atom for each subgraph. * \param ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled * \param GlobalAtomCount number of atoms in the complete molecule * \return true - success, false - failure (ListOfLocalAtoms != NULL) */ bool molecule::FillListOfLocalAtoms(ListOfLocalAtoms_t &ListOfLocalAtoms, const int GlobalAtomCount) { bool status = true; if (ListOfLocalAtoms.empty()) { // allocate and fill list of this fragment/subgraph status = status && CreateFatherLookupTable(ListOfLocalAtoms, GlobalAtomCount); } else return false; return status; } /** Creates a lookup table for true father's Atom::Nr -> atom ptr. * \param *start begin of list (STL iterator, i.e. first item) * \paran *end end of list (STL iterator, i.e. one past last item) * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start) * \param count optional predetermined size for table (otherwise we set the count to highest true father id) * \return true - success, false - failure */ bool molecule::CreateFatherLookupTable(ListOfLocalAtoms_t &LookupTable, int count) { bool status = true; int AtomNo; if (!LookupTable.empty()) { ELOG(1, "Pointer for Lookup table is not empty! Aborting ..."); return false; } // count them if (count == 0) { for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::Nr -> atom) used as a marker table lateron count = (count < (*iter)->GetTrueFather()->getNr()) ? (*iter)->GetTrueFather()->getNr() : count; } } if (count <= 0) { ELOG(1, "Count of lookup list is 0 or less."); return false; } // allocate and fill for (int i=0;i<=count;i++) LookupTable[i] = NULL; for (molecule::iterator iter = begin(); iter != end(); ++iter) { AtomNo = (*iter)->GetTrueFather()->getNr(); if ((AtomNo >= 0) && (AtomNo <= count)) { LOG(3, "DEBUG: Setting LookupTable[" << AtomNo << "] to " << *(*iter)); LookupTable[AtomNo] = (*iter); } else { ELOG(1, "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << "]."); status = false; break; } } return status; }; /** Corrects the nuclei position if the fragment was created over the cell borders. * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction. * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically * and re-add the bond. Looping on the distance check. * \param *out ofstream for debugging messages */ bool molecule::ScanForPeriodicCorrection() { bond::ptr Binder; //bond::ptr OtherBinder = NULL; atom *Walker = NULL; atom *OtherWalker = NULL; RealSpaceMatrix matrix = World::getInstance().getDomain().getM(); enum GraphEdge::Shading *ColorList = NULL; double tmp; //bool LastBond = true; // only needed to due list construct Vector Translationvector; //std::deque *CompStack = NULL; std::deque *AtomStack = new std::deque; // (getAtomCount()); bool flag = true; BondGraph *BG = World::getInstance().getBondGraph(); LOG(2, "Begin of ScanForPeriodicCorrection."); ColorList = new enum GraphEdge::Shading[getAtomCount()]; for (int i=0;igetListOfBonds(); for(BondList::const_iterator BondRunner = ListOfBonds.begin(); (!flag) && (BondRunner != ListOfBonds.end()); ++BondRunner) { Binder = (*BondRunner); for (int i=NDIM;i--;) { tmp = fabs(Binder->leftatom->at(i) - Binder->rightatom->at(i)); //LOG(3, "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "."); const range MinMaxDistance( BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom)); if (!MinMaxDistance.isInRange(tmp)) { LOG(2, "Correcting at bond " << *Binder << "."); flag = true; break; } } } } //if (flag) { if (0) { // create translation vector from their periodically modified distance for (int i=NDIM;i--;) { tmp = Binder->leftatom->at(i) - Binder->rightatom->at(i); const range MinMaxDistance( BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom)); if (fabs(tmp) > MinMaxDistance.last) // check against Min is not useful for components Translationvector[i] = (tmp < 0) ? +1. : -1.; } Translationvector *= matrix; LOG(3, "INFO: Translation vector is " << Translationvector << "."); // apply to all atoms of first component via BFS for (int i=getAtomCount();i--;) ColorList[i] = GraphEdge::white; AtomStack->push_front(Binder->leftatom); while (!AtomStack->empty()) { Walker = AtomStack->front(); AtomStack->pop_front(); //LOG(3, "INFO: Current Walker is: " << *Walker << "."); ColorList[Walker->getNr()] = GraphEdge::black; // mark as explored *Walker += Translationvector; // translate const BondList& ListOfBonds = Walker->getListOfBonds(); for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) { if ((*Runner) != Binder) { OtherWalker = (*Runner)->GetOtherAtom(Walker); if (ColorList[OtherWalker->getNr()] == GraphEdge::white) { AtomStack->push_front(OtherWalker); // push if yet unexplored } } } } // // re-add bond // if (OtherBinder == NULL) { // is the only bond? // //Do nothing // } else { // if (!LastBond) { // link(Binder, OtherBinder); // no more implemented bond::previous ... // } else { // link(OtherBinder, Binder); // no more implemented bond::previous ... // } // } } else { LOG(3, "No corrections for this fragment."); } //delete(CompStack); } // free allocated space from ReturnFullMatrixforSymmetric() delete(AtomStack); delete[](ColorList); LOG(2, "End of ScanForPeriodicCorrection."); return flag; };