| [bcf653] | 1 | /*
 | 
|---|
 | 2 |  * Project: MoleCuilder
 | 
|---|
 | 3 |  * Description: creates and alters molecular systems
 | 
|---|
 | 4 |  * Copyright (C)  2010 University of Bonn. All rights reserved.
 | 
|---|
 | 5 |  * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 | 
|---|
 | 6 |  */
 | 
|---|
 | 7 | 
 | 
|---|
| [cee0b57] | 8 | /*
 | 
|---|
 | 9 |  * molecule_graph.cpp
 | 
|---|
 | 10 |  *
 | 
|---|
 | 11 |  *  Created on: Oct 5, 2009
 | 
|---|
 | 12 |  *      Author: heber
 | 
|---|
 | 13 |  */
 | 
|---|
 | 14 | 
 | 
|---|
| [bf3817] | 15 | // include config.h
 | 
|---|
| [aafd77] | 16 | #ifdef HAVE_CONFIG_H
 | 
|---|
 | 17 | #include <config.h>
 | 
|---|
 | 18 | #endif
 | 
|---|
 | 19 | 
 | 
|---|
| [ad011c] | 20 | #include "CodePatterns/MemDebug.hpp"
 | 
|---|
| [112b09] | 21 | 
 | 
|---|
| [a564be] | 22 | #include <stack>
 | 
|---|
 | 23 | 
 | 
|---|
| [f66195] | 24 | #include "atom.hpp"
 | 
|---|
| [129204] | 25 | #include "Bond/bond.hpp"
 | 
|---|
| [34c43a] | 26 | #include "Box.hpp"
 | 
|---|
 | 27 | #include "CodePatterns/Assert.hpp"
 | 
|---|
 | 28 | #include "CodePatterns/Info.hpp"
 | 
|---|
 | 29 | #include "CodePatterns/Log.hpp"
 | 
|---|
 | 30 | #include "CodePatterns/Verbose.hpp"
 | 
|---|
| [cee0b57] | 31 | #include "config.hpp"
 | 
|---|
| [49c059] | 32 | #include "Graph/DepthFirstSearchAnalysis.hpp"
 | 
|---|
| [f66195] | 33 | #include "element.hpp"
 | 
|---|
| [129204] | 34 | #include "Graph/BondGraph.hpp"
 | 
|---|
| [1d5afa5] | 35 | #include "Helpers/defs.hpp"
 | 
|---|
 | 36 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
 | 
|---|
| [b8b75d] | 37 | #include "linkedcell.hpp"
 | 
|---|
| [cee0b57] | 38 | #include "molecule.hpp"
 | 
|---|
| [34c43a] | 39 | #include "PointCloudAdaptor.hpp"
 | 
|---|
| [b34306] | 40 | #include "World.hpp"
 | 
|---|
| [9d83b6] | 41 | #include "WorldTime.hpp"
 | 
|---|
| [1d5afa5] | 42 | 
 | 
|---|
| [cee0b57] | 43 | 
 | 
|---|
| [99752a] | 44 | /** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
 | 
|---|
 | 45 |  * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
 | 
|---|
 | 46 |  * \param *reference reference molecule with the bond structure to be copied
 | 
|---|
 | 47 |  * \param **&ListOfLocalAtoms Lookup table for this subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
 | 
|---|
 | 48 |  * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 | 
|---|
 | 49 |  * \return true - success, false - failure
 | 
|---|
 | 50 |  */
 | 
|---|
 | 51 | bool molecule::FillBondStructureFromReference(const molecule * const reference, atom **&ListOfLocalAtoms, bool FreeList)
 | 
|---|
 | 52 | {
 | 
|---|
 | 53 |   atom *OtherWalker = NULL;
 | 
|---|
 | 54 |   atom *Father = NULL;
 | 
|---|
 | 55 |   bool status = true;
 | 
|---|
 | 56 |   int AtomNo;
 | 
|---|
 | 57 | 
 | 
|---|
 | 58 |   DoLog(1) && (Log() << Verbose(1) << "Begin of FillBondStructureFromReference." << endl);
 | 
|---|
 | 59 |   // fill ListOfLocalAtoms if NULL was given
 | 
|---|
 | 60 |   if (!FillListOfLocalAtoms(ListOfLocalAtoms, reference->getAtomCount())) {
 | 
|---|
 | 61 |     DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl);
 | 
|---|
 | 62 |     return false;
 | 
|---|
 | 63 |   }
 | 
|---|
 | 64 | 
 | 
|---|
 | 65 |   if (status) {
 | 
|---|
 | 66 |     DoLog(1) && (Log() << Verbose(1) << "Creating adjacency list for molecule " << getName() << "." << endl);
 | 
|---|
 | 67 |     // remove every bond from the list
 | 
|---|
 | 68 |     for_each(begin(), end(),
 | 
|---|
 | 69 |         boost::bind(&BondedParticle::ClearBondsAtStep,_1,WorldTime::getTime()));
 | 
|---|
 | 70 | 
 | 
|---|
 | 71 | 
 | 
|---|
 | 72 |     for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 | 
|---|
 | 73 |       Father = (*iter)->GetTrueFather();
 | 
|---|
 | 74 |       AtomNo = Father->getNr(); // global id of the current walker
 | 
|---|
 | 75 |       const BondList& ListOfBonds = Father->getListOfBonds();
 | 
|---|
 | 76 |       for (BondList::const_iterator Runner = ListOfBonds.begin();
 | 
|---|
 | 77 |           Runner != ListOfBonds.end();
 | 
|---|
 | 78 |           ++Runner) {
 | 
|---|
 | 79 |         OtherWalker = ListOfLocalAtoms[(*Runner)->GetOtherAtom((*iter)->GetTrueFather())->getNr()]; // local copy of current bond partner of walker
 | 
|---|
 | 80 |         if (OtherWalker != NULL) {
 | 
|---|
 | 81 |           if (OtherWalker->getNr() > (*iter)->getNr())
 | 
|---|
 | 82 |             AddBond((*iter), OtherWalker, (*Runner)->BondDegree);
 | 
|---|
 | 83 |         } else {
 | 
|---|
 | 84 |           DoLog(1) && (Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << (*Runner)->GetOtherAtom((*iter)->GetTrueFather())->getNr() << "] is NULL!" << endl);
 | 
|---|
 | 85 |           status = false;
 | 
|---|
 | 86 |         }
 | 
|---|
 | 87 |       }
 | 
|---|
 | 88 |     }
 | 
|---|
 | 89 |   }
 | 
|---|
 | 90 | 
 | 
|---|
 | 91 |   if ((FreeList) && (ListOfLocalAtoms != NULL)) {
 | 
|---|
 | 92 |     // free the index lookup list
 | 
|---|
 | 93 |     delete[](ListOfLocalAtoms);
 | 
|---|
 | 94 |   }
 | 
|---|
 | 95 |   DoLog(1) && (Log() << Verbose(1) << "End of FillBondStructureFromReference." << endl);
 | 
|---|
 | 96 |   return status;
 | 
|---|
 | 97 | };
 | 
|---|
 | 98 | 
 | 
|---|
| [e08c46] | 99 | /** Checks for presence of bonds within atom list.
 | 
|---|
 | 100 |  * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
 | 
|---|
 | 101 |  * \return true - bonds present, false - no bonds
 | 
|---|
 | 102 |  */
 | 
|---|
| [e4afb4] | 103 | bool molecule::hasBondStructure() const
 | 
|---|
| [e08c46] | 104 | {
 | 
|---|
| [9d83b6] | 105 |   for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
 | 
|---|
| [5e2f80] | 106 |     //LOG(0, "molecule::hasBondStructure() - checking bond list of atom " << (*AtomRunner)->getId() << ".");
 | 
|---|
| [9d83b6] | 107 |     const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
 | 
|---|
 | 108 |     if (!ListOfBonds.empty())
 | 
|---|
| [e08c46] | 109 |       return true;
 | 
|---|
| [9d83b6] | 110 |   }
 | 
|---|
| [e08c46] | 111 |   return false;
 | 
|---|
 | 112 | }
 | 
|---|
 | 113 | 
 | 
|---|
| [b8b75d] | 114 | /** Prints a list of all bonds to \a *out.
 | 
|---|
 | 115 |  */
 | 
|---|
| [e138de] | 116 | void molecule::OutputBondsList() const
 | 
|---|
| [b8b75d] | 117 | {
 | 
|---|
| [a67d19] | 118 |   DoLog(1) && (Log() << Verbose(1) << endl << "From contents of bond chain list:");
 | 
|---|
| [9d83b6] | 119 |   for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) {
 | 
|---|
 | 120 |     const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
 | 
|---|
 | 121 |     for(BondList::const_iterator BondRunner = ListOfBonds.begin();
 | 
|---|
 | 122 |         BondRunner != ListOfBonds.end();
 | 
|---|
 | 123 |         ++BondRunner)
 | 
|---|
| [e08c46] | 124 |       if ((*BondRunner)->leftatom == *AtomRunner) {
 | 
|---|
 | 125 |         DoLog(0) && (Log() << Verbose(0) << *(*BondRunner) << "\t" << endl);
 | 
|---|
 | 126 |       }
 | 
|---|
| [9d83b6] | 127 |   }
 | 
|---|
| [a67d19] | 128 |   DoLog(0) && (Log() << Verbose(0) << endl);
 | 
|---|
| [9eefda] | 129 | }
 | 
|---|
 | 130 | ;
 | 
|---|
| [cee0b57] | 131 | 
 | 
|---|
 | 132 | 
 | 
|---|
 | 133 | /** Storing the bond structure of a molecule to file.
 | 
|---|
| [5309ba] | 134 |  * Simply stores Atom::Nr and then the Atom::Nr of all bond partners per line.
 | 
|---|
| [35b698] | 135 |  * \param &filename name of file
 | 
|---|
 | 136 |  * \param path path to file, defaults to empty
 | 
|---|
| [cee0b57] | 137 |  * \return true - file written successfully, false - writing failed
 | 
|---|
 | 138 |  */
 | 
|---|
| [e4afb4] | 139 | bool molecule::StoreAdjacencyToFile(std::string filename, std::string path)
 | 
|---|
| [cee0b57] | 140 | {
 | 
|---|
 | 141 |   ofstream AdjacencyFile;
 | 
|---|
| [35b698] | 142 |   string line;
 | 
|---|
| [cee0b57] | 143 |   bool status = true;
 | 
|---|
 | 144 | 
 | 
|---|
| [35b698] | 145 |   if (path != "")
 | 
|---|
 | 146 |     line = path + "/" + filename;
 | 
|---|
| [8ab0407] | 147 |   else
 | 
|---|
| [35b698] | 148 |     line = filename;
 | 
|---|
 | 149 |   AdjacencyFile.open(line.c_str(), ios::out);
 | 
|---|
| [acf800] | 150 |   DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
 | 
|---|
| [35b698] | 151 |   if (AdjacencyFile.good()) {
 | 
|---|
| [1f1b23] | 152 |     AdjacencyFile << "m\tn" << endl;
 | 
|---|
| [00ef5c] | 153 |     for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputAdjacency),&AdjacencyFile));
 | 
|---|
| [cee0b57] | 154 |     AdjacencyFile.close();
 | 
|---|
| [acf800] | 155 |     DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
 | 
|---|
| [cee0b57] | 156 |   } else {
 | 
|---|
| [35b698] | 157 |     DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
 | 
|---|
| [cee0b57] | 158 |     status = false;
 | 
|---|
 | 159 |   }
 | 
|---|
 | 160 | 
 | 
|---|
 | 161 |   return status;
 | 
|---|
| [9eefda] | 162 | }
 | 
|---|
 | 163 | ;
 | 
|---|
| [cee0b57] | 164 | 
 | 
|---|
| [1f1b23] | 165 | /** Storing the bond structure of a molecule to file.
 | 
|---|
| [5309ba] | 166 |  * Simply stores Atom::Nr and then the Atom::Nr of all bond partners, one per line.
 | 
|---|
| [35b698] | 167 |  * \param &filename name of file
 | 
|---|
 | 168 |  * \param path path to file, defaults to empty
 | 
|---|
| [1f1b23] | 169 |  * \return true - file written successfully, false - writing failed
 | 
|---|
 | 170 |  */
 | 
|---|
| [e4afb4] | 171 | bool molecule::StoreBondsToFile(std::string filename, std::string path)
 | 
|---|
| [1f1b23] | 172 | {
 | 
|---|
 | 173 |   ofstream BondFile;
 | 
|---|
| [35b698] | 174 |   string line;
 | 
|---|
| [1f1b23] | 175 |   bool status = true;
 | 
|---|
 | 176 | 
 | 
|---|
| [35b698] | 177 |   if (path != "")
 | 
|---|
 | 178 |     line = path + "/" + filename;
 | 
|---|
| [8ab0407] | 179 |   else
 | 
|---|
| [35b698] | 180 |     line = filename;
 | 
|---|
 | 181 |   BondFile.open(line.c_str(), ios::out);
 | 
|---|
| [acf800] | 182 |   DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
 | 
|---|
| [35b698] | 183 |   if (BondFile.good()) {
 | 
|---|
| [1f1b23] | 184 |     BondFile << "m\tn" << endl;
 | 
|---|
| [00ef5c] | 185 |     for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile));
 | 
|---|
| [1f1b23] | 186 |     BondFile.close();
 | 
|---|
| [acf800] | 187 |     DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
 | 
|---|
| [1f1b23] | 188 |   } else {
 | 
|---|
| [35b698] | 189 |     DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
 | 
|---|
| [1f1b23] | 190 |     status = false;
 | 
|---|
 | 191 |   }
 | 
|---|
 | 192 | 
 | 
|---|
 | 193 |   return status;
 | 
|---|
 | 194 | }
 | 
|---|
 | 195 | ;
 | 
|---|
 | 196 | 
 | 
|---|
| [266237] | 197 | /** Adds a bond as a copy to a given one
 | 
|---|
 | 198 |  * \param *left leftatom of new bond
 | 
|---|
 | 199 |  * \param *right rightatom of new bond
 | 
|---|
 | 200 |  * \param *CopyBond rest of fields in bond are copied from this
 | 
|---|
 | 201 |  * \return pointer to new bond
 | 
|---|
 | 202 |  */
 | 
|---|
 | 203 | bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
 | 
|---|
 | 204 | {
 | 
|---|
 | 205 |   bond *Binder = AddBond(left, right, CopyBond->BondDegree);
 | 
|---|
 | 206 |   Binder->Cyclic = CopyBond->Cyclic;
 | 
|---|
 | 207 |   Binder->Type = CopyBond->Type;
 | 
|---|
 | 208 |   return Binder;
 | 
|---|
| [9eefda] | 209 | }
 | 
|---|
 | 210 | ;
 | 
|---|
| [266237] | 211 | 
 | 
|---|
| [49c059] | 212 | /** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
 | 
|---|
| [c6123b] | 213 |  * \param **&ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 | 
|---|
 | 214 |  * \param GlobalAtomCount number of atoms in the complete molecule
 | 
|---|
 | 215 |  * \return true - success, false - failure (ListOfLocalAtoms != NULL)
 | 
|---|
 | 216 |  */
 | 
|---|
 | 217 | bool molecule::FillListOfLocalAtoms(atom **&ListOfLocalAtoms, const int GlobalAtomCount)
 | 
|---|
 | 218 | {
 | 
|---|
 | 219 |   bool status = true;
 | 
|---|
 | 220 | 
 | 
|---|
 | 221 |   if (ListOfLocalAtoms == NULL) { // allocate and fill list of this fragment/subgraph
 | 
|---|
 | 222 |     status = status && CreateFatherLookupTable(ListOfLocalAtoms, GlobalAtomCount);
 | 
|---|
 | 223 |   } else
 | 
|---|
 | 224 |     return false;
 | 
|---|
 | 225 | 
 | 
|---|
 | 226 |   return status;
 | 
|---|
 | 227 | }
 | 
|---|
 | 228 | 
 | 
|---|