[bcf653] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
[0aa122] | 4 | * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
|
---|
[5aaa43] | 5 | * Copyright (C) 2013 Frederik Heber. All rights reserved.
|
---|
[94d5ac6] | 6 | *
|
---|
| 7 | *
|
---|
| 8 | * This file is part of MoleCuilder.
|
---|
| 9 | *
|
---|
| 10 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 11 | * it under the terms of the GNU General Public License as published by
|
---|
| 12 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 13 | * (at your option) any later version.
|
---|
| 14 | *
|
---|
| 15 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 16 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 17 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 18 | * GNU General Public License for more details.
|
---|
| 19 | *
|
---|
| 20 | * You should have received a copy of the GNU General Public License
|
---|
| 21 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
[bcf653] | 22 | */
|
---|
| 23 |
|
---|
[cee0b57] | 24 | /*
|
---|
| 25 | * molecule_graph.cpp
|
---|
| 26 | *
|
---|
| 27 | * Created on: Oct 5, 2009
|
---|
| 28 | * Author: heber
|
---|
| 29 | */
|
---|
| 30 |
|
---|
[bf3817] | 31 | // include config.h
|
---|
[aafd77] | 32 | #ifdef HAVE_CONFIG_H
|
---|
| 33 | #include <config.h>
|
---|
| 34 | #endif
|
---|
| 35 |
|
---|
[9eb71b3] | 36 | //#include "CodePatterns/MemDebug.hpp"
|
---|
[112b09] | 37 |
|
---|
[a564be] | 38 | #include <stack>
|
---|
| 39 |
|
---|
[6f0841] | 40 | #include "Atom/atom.hpp"
|
---|
[129204] | 41 | #include "Bond/bond.hpp"
|
---|
[34c43a] | 42 | #include "Box.hpp"
|
---|
| 43 | #include "CodePatterns/Assert.hpp"
|
---|
| 44 | #include "CodePatterns/Info.hpp"
|
---|
| 45 | #include "CodePatterns/Log.hpp"
|
---|
| 46 | #include "CodePatterns/Verbose.hpp"
|
---|
[cee0b57] | 47 | #include "config.hpp"
|
---|
[49c059] | 48 | #include "Graph/DepthFirstSearchAnalysis.hpp"
|
---|
[3bdb6d] | 49 | #include "Element/element.hpp"
|
---|
[129204] | 50 | #include "Graph/BondGraph.hpp"
|
---|
[6d551c] | 51 | #include "Graph/ListOfLocalAtoms.hpp"
|
---|
[1d5afa5] | 52 | #include "Helpers/defs.hpp"
|
---|
[246e13] | 53 | #include "Helpers/helpers.hpp"
|
---|
[1d5afa5] | 54 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
[53c7fc] | 55 | #include "LinkedCell/linkedcell.hpp"
|
---|
| 56 | #include "LinkedCell/PointCloudAdaptor.hpp"
|
---|
[cee0b57] | 57 | #include "molecule.hpp"
|
---|
[b34306] | 58 | #include "World.hpp"
|
---|
[9d83b6] | 59 | #include "WorldTime.hpp"
|
---|
[1d5afa5] | 60 |
|
---|
[cee0b57] | 61 |
|
---|
[99752a] | 62 | /** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
|
---|
| 63 | * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
|
---|
| 64 | * \param *reference reference molecule with the bond structure to be copied
|
---|
[6d551c] | 65 | * \param ListOfLocalAtoms Lookup table for this subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
|
---|
| 66 | * \param FreeList true - ListOfLocalAtoms is free'd before return, false - it is not
|
---|
[99752a] | 67 | * \return true - success, false - failure
|
---|
| 68 | */
|
---|
[6d551c] | 69 | bool molecule::FillBondStructureFromReference(const molecule * const reference, ListOfLocalAtoms_t &ListOfLocalAtoms, bool FreeList)
|
---|
[99752a] | 70 | {
|
---|
| 71 | bool status = true;
|
---|
| 72 |
|
---|
[47d041] | 73 | LOG(1, "Begin of FillBondStructureFromReference.");
|
---|
[99752a] | 74 | // fill ListOfLocalAtoms if NULL was given
|
---|
| 75 | if (!FillListOfLocalAtoms(ListOfLocalAtoms, reference->getAtomCount())) {
|
---|
[47d041] | 76 | LOG(1, "Filling of ListOfLocalAtoms failed.");
|
---|
[99752a] | 77 | return false;
|
---|
| 78 | }
|
---|
| 79 |
|
---|
| 80 | if (status) {
|
---|
[47d041] | 81 | LOG(1, "Creating adjacency list for molecule " << getName() << ".");
|
---|
[99752a] | 82 | // remove every bond from the list
|
---|
| 83 | for_each(begin(), end(),
|
---|
| 84 | boost::bind(&BondedParticle::ClearBondsAtStep,_1,WorldTime::getTime()));
|
---|
| 85 |
|
---|
| 86 |
|
---|
[59fff1] | 87 | for(molecule::iterator iter = begin(); iter != end(); ++iter) {
|
---|
| 88 | const atom * const Father = (*iter)->GetTrueFather();
|
---|
[c9cafa] | 89 | //const int AtomNo = Father->getNr(); // global id of the current walker
|
---|
[99752a] | 90 | const BondList& ListOfBonds = Father->getListOfBonds();
|
---|
| 91 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 92 | Runner != ListOfBonds.end();
|
---|
| 93 | ++Runner) {
|
---|
[59fff1] | 94 | atom * const OtherAtom = (*Runner)->GetOtherAtom((*iter)->GetTrueFather());
|
---|
[6d551c] | 95 | const ListOfLocalAtoms_t::const_iterator localiter = ListOfLocalAtoms.find(OtherAtom->getNr());
|
---|
| 96 | ASSERT( localiter != ListOfLocalAtoms.end(),
|
---|
| 97 | "molecule::FillBondStructureFromReference() - could not find id"
|
---|
| 98 | +toString(OtherAtom->getNr())+" in ListOfLocalAtoms.");
|
---|
| 99 | atom * const OtherWalker = localiter->second; // local copy of current bond partner of walker
|
---|
[99752a] | 100 | if (OtherWalker != NULL) {
|
---|
| 101 | if (OtherWalker->getNr() > (*iter)->getNr())
|
---|
[1f693d] | 102 | AddBond((*iter), OtherWalker, (*Runner)->getDegree());
|
---|
[99752a] | 103 | } else {
|
---|
[59fff1] | 104 | LOG(1, "OtherWalker = ListOfLocalAtoms[" << OtherAtom->getNr() << "] is NULL!");
|
---|
[99752a] | 105 | status = false;
|
---|
| 106 | }
|
---|
| 107 | }
|
---|
| 108 | }
|
---|
| 109 | }
|
---|
| 110 |
|
---|
[6d551c] | 111 | if ((FreeList) && (!ListOfLocalAtoms.empty())) {
|
---|
[99752a] | 112 | // free the index lookup list
|
---|
[6d551c] | 113 | ListOfLocalAtoms.clear();
|
---|
[99752a] | 114 | }
|
---|
[47d041] | 115 | LOG(1, "End of FillBondStructureFromReference.");
|
---|
[99752a] | 116 | return status;
|
---|
| 117 | };
|
---|
| 118 |
|
---|
[e08c46] | 119 | /** Checks for presence of bonds within atom list.
|
---|
| 120 | * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
|
---|
| 121 | * \return true - bonds present, false - no bonds
|
---|
| 122 | */
|
---|
[e4afb4] | 123 | bool molecule::hasBondStructure() const
|
---|
[e08c46] | 124 | {
|
---|
[9d83b6] | 125 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
[5e2f80] | 126 | //LOG(0, "molecule::hasBondStructure() - checking bond list of atom " << (*AtomRunner)->getId() << ".");
|
---|
[9d83b6] | 127 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 128 | if (!ListOfBonds.empty())
|
---|
[e08c46] | 129 | return true;
|
---|
[9d83b6] | 130 | }
|
---|
[e08c46] | 131 | return false;
|
---|
| 132 | }
|
---|
| 133 |
|
---|
[b8b75d] | 134 | /** Prints a list of all bonds to \a *out.
|
---|
| 135 | */
|
---|
[e138de] | 136 | void molecule::OutputBondsList() const
|
---|
[b8b75d] | 137 | {
|
---|
[47d041] | 138 | if (DoLog(1)) {
|
---|
| 139 | std::stringstream output;
|
---|
| 140 | output << std::endl << "From contents of bond chain list:";
|
---|
| 141 | for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) {
|
---|
| 142 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 143 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 144 | BondRunner != ListOfBonds.end();
|
---|
| 145 | ++BondRunner)
|
---|
| 146 | if ((*BondRunner)->leftatom == *AtomRunner) {
|
---|
| 147 | output << *(*BondRunner) << "\t";
|
---|
| 148 | }
|
---|
| 149 | }
|
---|
| 150 | LOG(1, output.str());
|
---|
[9d83b6] | 151 | }
|
---|
[9eefda] | 152 | }
|
---|
[cee0b57] | 153 |
|
---|
| 154 |
|
---|
[1f1b23] | 155 | /** Storing the bond structure of a molecule to file.
|
---|
[5309ba] | 156 | * Simply stores Atom::Nr and then the Atom::Nr of all bond partners, one per line.
|
---|
[35b698] | 157 | * \param &filename name of file
|
---|
| 158 | * \param path path to file, defaults to empty
|
---|
[1f1b23] | 159 | * \return true - file written successfully, false - writing failed
|
---|
| 160 | */
|
---|
[e4afb4] | 161 | bool molecule::StoreBondsToFile(std::string filename, std::string path)
|
---|
[1f1b23] | 162 | {
|
---|
| 163 | ofstream BondFile;
|
---|
[35b698] | 164 | string line;
|
---|
[1f1b23] | 165 | bool status = true;
|
---|
| 166 |
|
---|
[35b698] | 167 | if (path != "")
|
---|
| 168 | line = path + "/" + filename;
|
---|
[8ab0407] | 169 | else
|
---|
[35b698] | 170 | line = filename;
|
---|
| 171 | BondFile.open(line.c_str(), ios::out);
|
---|
[47d041] | 172 | LOG(1, "Saving adjacency list ... ");
|
---|
[35b698] | 173 | if (BondFile.good()) {
|
---|
[1f1b23] | 174 | BondFile << "m\tn" << endl;
|
---|
[30c753] | 175 | for_each(begin(),end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile));
|
---|
[1f1b23] | 176 | BondFile.close();
|
---|
[47d041] | 177 | LOG(1, "\t... done.");
|
---|
[1f1b23] | 178 | } else {
|
---|
[47d041] | 179 | LOG(1, "\t... failed to open file " << line << ".");
|
---|
[1f1b23] | 180 | status = false;
|
---|
| 181 | }
|
---|
| 182 |
|
---|
| 183 | return status;
|
---|
| 184 | }
|
---|
| 185 | ;
|
---|
| 186 |
|
---|
[266237] | 187 | /** Adds a bond as a copy to a given one
|
---|
| 188 | * \param *left leftatom of new bond
|
---|
| 189 | * \param *right rightatom of new bond
|
---|
| 190 | * \param *CopyBond rest of fields in bond are copied from this
|
---|
| 191 | * \return pointer to new bond
|
---|
| 192 | */
|
---|
[88c8ec] | 193 | bond::ptr molecule::CopyBond(atom *left, atom *right, bond::ptr CopyBond)
|
---|
[266237] | 194 | {
|
---|
[1f693d] | 195 | bond::ptr Binder = AddBond(left, right, CopyBond->getDegree());
|
---|
[266237] | 196 | Binder->Cyclic = CopyBond->Cyclic;
|
---|
| 197 | Binder->Type = CopyBond->Type;
|
---|
| 198 | return Binder;
|
---|
[9eefda] | 199 | }
|
---|
| 200 | ;
|
---|
[266237] | 201 |
|
---|
[49c059] | 202 | /** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
|
---|
[6d551c] | 203 | * \param ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
|
---|
[c6123b] | 204 | * \param GlobalAtomCount number of atoms in the complete molecule
|
---|
| 205 | * \return true - success, false - failure (ListOfLocalAtoms != NULL)
|
---|
| 206 | */
|
---|
[6d551c] | 207 | bool molecule::FillListOfLocalAtoms(ListOfLocalAtoms_t &ListOfLocalAtoms, const int GlobalAtomCount)
|
---|
[c6123b] | 208 | {
|
---|
| 209 | bool status = true;
|
---|
| 210 |
|
---|
[6d551c] | 211 | if (ListOfLocalAtoms.empty()) { // allocate and fill list of this fragment/subgraph
|
---|
[c6123b] | 212 | status = status && CreateFatherLookupTable(ListOfLocalAtoms, GlobalAtomCount);
|
---|
| 213 | } else
|
---|
| 214 | return false;
|
---|
| 215 |
|
---|
| 216 | return status;
|
---|
| 217 | }
|
---|
| 218 |
|
---|
[246e13] | 219 |
|
---|
| 220 | /** Creates a lookup table for true father's Atom::Nr -> atom ptr.
|
---|
| 221 | * \param *start begin of list (STL iterator, i.e. first item)
|
---|
| 222 | * \paran *end end of list (STL iterator, i.e. one past last item)
|
---|
| 223 | * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
|
---|
| 224 | * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
|
---|
| 225 | * \return true - success, false - failure
|
---|
| 226 | */
|
---|
[6d551c] | 227 | bool molecule::CreateFatherLookupTable(ListOfLocalAtoms_t &LookupTable, int count)
|
---|
[246e13] | 228 | {
|
---|
| 229 | bool status = true;
|
---|
| 230 | int AtomNo;
|
---|
| 231 |
|
---|
[6d551c] | 232 | if (!LookupTable.empty()) {
|
---|
| 233 | ELOG(1, "Pointer for Lookup table is not empty! Aborting ...");
|
---|
[246e13] | 234 | return false;
|
---|
| 235 | }
|
---|
| 236 |
|
---|
| 237 | // count them
|
---|
| 238 | if (count == 0) {
|
---|
| 239 | for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::Nr -> atom) used as a marker table lateron
|
---|
| 240 | count = (count < (*iter)->GetTrueFather()->getNr()) ? (*iter)->GetTrueFather()->getNr() : count;
|
---|
| 241 | }
|
---|
| 242 | }
|
---|
| 243 | if (count <= 0) {
|
---|
[47d041] | 244 | ELOG(1, "Count of lookup list is 0 or less.");
|
---|
[246e13] | 245 | return false;
|
---|
| 246 | }
|
---|
| 247 |
|
---|
| 248 | // allocate and fill
|
---|
[6d551c] | 249 | for (int i=0;i<=count;i++)
|
---|
| 250 | LookupTable[i] = NULL;
|
---|
| 251 | for (molecule::iterator iter = begin(); iter != end(); ++iter) {
|
---|
| 252 | AtomNo = (*iter)->GetTrueFather()->getNr();
|
---|
| 253 | if ((AtomNo >= 0) && (AtomNo <= count)) {
|
---|
| 254 | LOG(3, "DEBUG: Setting LookupTable[" << AtomNo << "] to " << *(*iter));
|
---|
| 255 | LookupTable[AtomNo] = (*iter);
|
---|
| 256 | } else {
|
---|
| 257 | ELOG(1, "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << "].");
|
---|
| 258 | status = false;
|
---|
| 259 | break;
|
---|
[246e13] | 260 | }
|
---|
| 261 | }
|
---|
| 262 |
|
---|
| 263 | return status;
|
---|
| 264 | };
|
---|
| 265 |
|
---|
| 266 |
|
---|
| 267 |
|
---|
| 268 | /** Corrects the nuclei position if the fragment was created over the cell borders.
|
---|
| 269 | * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
|
---|
| 270 | * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
|
---|
| 271 | * and re-add the bond. Looping on the distance check.
|
---|
| 272 | * \param *out ofstream for debugging messages
|
---|
| 273 | */
|
---|
| 274 | bool molecule::ScanForPeriodicCorrection()
|
---|
| 275 | {
|
---|
[7d82a5] | 276 | bond::ptr Binder;
|
---|
[88c8ec] | 277 | //bond::ptr OtherBinder = NULL;
|
---|
[246e13] | 278 | atom *Walker = NULL;
|
---|
| 279 | atom *OtherWalker = NULL;
|
---|
| 280 | RealSpaceMatrix matrix = World::getInstance().getDomain().getM();
|
---|
| 281 | enum GraphEdge::Shading *ColorList = NULL;
|
---|
| 282 | double tmp;
|
---|
| 283 | //bool LastBond = true; // only needed to due list construct
|
---|
| 284 | Vector Translationvector;
|
---|
| 285 | //std::deque<atom *> *CompStack = NULL;
|
---|
| 286 | std::deque<atom *> *AtomStack = new std::deque<atom *>; // (getAtomCount());
|
---|
| 287 | bool flag = true;
|
---|
| 288 | BondGraph *BG = World::getInstance().getBondGraph();
|
---|
| 289 |
|
---|
[47d041] | 290 | LOG(2, "Begin of ScanForPeriodicCorrection.");
|
---|
[246e13] | 291 |
|
---|
| 292 | ColorList = new enum GraphEdge::Shading[getAtomCount()];
|
---|
| 293 | for (int i=0;i<getAtomCount();i++)
|
---|
| 294 | ColorList[i] = (enum GraphEdge::Shading)0;
|
---|
| 295 | if (flag) {
|
---|
| 296 | // remove bonds that are beyond bonddistance
|
---|
| 297 | Translationvector.Zero();
|
---|
| 298 | // scan all bonds
|
---|
| 299 | flag = false;
|
---|
| 300 | for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner) {
|
---|
| 301 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 302 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 303 | (!flag) && (BondRunner != ListOfBonds.end());
|
---|
| 304 | ++BondRunner) {
|
---|
| 305 | Binder = (*BondRunner);
|
---|
| 306 | for (int i=NDIM;i--;) {
|
---|
| 307 | tmp = fabs(Binder->leftatom->at(i) - Binder->rightatom->at(i));
|
---|
[47d041] | 308 | //LOG(3, "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << ".");
|
---|
[246e13] | 309 | const range<double> MinMaxDistance(
|
---|
| 310 | BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
|
---|
| 311 | if (!MinMaxDistance.isInRange(tmp)) {
|
---|
[47d041] | 312 | LOG(2, "Correcting at bond " << *Binder << ".");
|
---|
[246e13] | 313 | flag = true;
|
---|
| 314 | break;
|
---|
| 315 | }
|
---|
| 316 | }
|
---|
| 317 | }
|
---|
| 318 | }
|
---|
| 319 | //if (flag) {
|
---|
| 320 | if (0) {
|
---|
| 321 | // create translation vector from their periodically modified distance
|
---|
| 322 | for (int i=NDIM;i--;) {
|
---|
| 323 | tmp = Binder->leftatom->at(i) - Binder->rightatom->at(i);
|
---|
| 324 | const range<double> MinMaxDistance(
|
---|
| 325 | BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
|
---|
| 326 | if (fabs(tmp) > MinMaxDistance.last) // check against Min is not useful for components
|
---|
| 327 | Translationvector[i] = (tmp < 0) ? +1. : -1.;
|
---|
| 328 | }
|
---|
| 329 | Translationvector *= matrix;
|
---|
[47d041] | 330 | LOG(3, "INFO: Translation vector is " << Translationvector << ".");
|
---|
[246e13] | 331 | // apply to all atoms of first component via BFS
|
---|
| 332 | for (int i=getAtomCount();i--;)
|
---|
| 333 | ColorList[i] = GraphEdge::white;
|
---|
| 334 | AtomStack->push_front(Binder->leftatom);
|
---|
| 335 | while (!AtomStack->empty()) {
|
---|
| 336 | Walker = AtomStack->front();
|
---|
| 337 | AtomStack->pop_front();
|
---|
[47d041] | 338 | //LOG(3, "INFO: Current Walker is: " << *Walker << ".");
|
---|
[246e13] | 339 | ColorList[Walker->getNr()] = GraphEdge::black; // mark as explored
|
---|
| 340 | *Walker += Translationvector; // translate
|
---|
| 341 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 342 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 343 | Runner != ListOfBonds.end();
|
---|
| 344 | ++Runner) {
|
---|
| 345 | if ((*Runner) != Binder) {
|
---|
| 346 | OtherWalker = (*Runner)->GetOtherAtom(Walker);
|
---|
| 347 | if (ColorList[OtherWalker->getNr()] == GraphEdge::white) {
|
---|
| 348 | AtomStack->push_front(OtherWalker); // push if yet unexplored
|
---|
| 349 | }
|
---|
| 350 | }
|
---|
| 351 | }
|
---|
| 352 | }
|
---|
| 353 | // // re-add bond
|
---|
| 354 | // if (OtherBinder == NULL) { // is the only bond?
|
---|
| 355 | // //Do nothing
|
---|
| 356 | // } else {
|
---|
| 357 | // if (!LastBond) {
|
---|
| 358 | // link(Binder, OtherBinder); // no more implemented bond::previous ...
|
---|
| 359 | // } else {
|
---|
| 360 | // link(OtherBinder, Binder); // no more implemented bond::previous ...
|
---|
| 361 | // }
|
---|
| 362 | // }
|
---|
| 363 | } else {
|
---|
[47d041] | 364 | LOG(3, "No corrections for this fragment.");
|
---|
[246e13] | 365 | }
|
---|
| 366 | //delete(CompStack);
|
---|
| 367 | }
|
---|
| 368 | // free allocated space from ReturnFullMatrixforSymmetric()
|
---|
| 369 | delete(AtomStack);
|
---|
| 370 | delete[](ColorList);
|
---|
[47d041] | 371 | LOG(2, "End of ScanForPeriodicCorrection.");
|
---|
[246e13] | 372 |
|
---|
| 373 | return flag;
|
---|
| 374 | };
|
---|