/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2012 University of Bonn. 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 . */ /* * ExportGraph_ToFiles.cpp * * Created on: 08.03.2012 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "ExportGraph_ToFiles.hpp" #include "CodePatterns/Info.hpp" #include "CodePatterns/Log.hpp" #include "Bond/bond.hpp" #include "Element/element.hpp" #include "Fragmentation/Graph.hpp" #include "Fragmentation/KeySet.hpp" #include "Fragmentation/SortIndex.hpp" #include "Graph/ListOfLocalAtoms.hpp" #include "molecule.hpp" #include "MoleculeListClass.hpp" #include "Parser/FormatParserStorage.hpp" #include "World.hpp" /** Constructor for class ExportGraph_ToFiles. * * @param _graph instance of Graph containing keyset of each fragment * @param saturation whether to saturate dangling bonds with hydrogen or not */ ExportGraph_ToFiles::ExportGraph_ToFiles(const Graph &_graph, const enum HydrogenSaturation _saturation) : ExportGraph(_graph), BondFragments(World::getPointer()), saturation(_saturation) {} /** Destructor of class ExportGraph_ToFiles. * * We free all created molecules again and also removed their copied atoms. */ ExportGraph_ToFiles::~ExportGraph_ToFiles() { // remove all create molecules again from the World including their atoms for (MoleculeList::iterator iter = BondFragments.ListOfMolecules.begin(); !BondFragments.ListOfMolecules.empty(); iter = BondFragments.ListOfMolecules.begin()) { // remove copied atoms and molecule again molecule *mol = *iter; mol->removeAtomsinMolecule(); World::getInstance().destroyMolecule(mol); BondFragments.ListOfMolecules.erase(iter); } } /** Actual implementation of the export to files function. */ void ExportGraph_ToFiles::operator()() { if (BondFragments.ListOfMolecules.size() == 0) prepareMolecule(); // ===== 9. Save fragments' configuration and keyset files et al to disk === LOG(1, "Writing " << BondFragments.ListOfMolecules.size() << " possible bond fragmentation configs"); bool write_status = true; for (std::vector::const_iterator iter = typelist.begin(); iter != typelist.end(); ++iter) { LOG(2, "INFO: Writing bond fragments for type " << (*iter) << "."); write_status = write_status && BondFragments.OutputConfigForListOfFragments( prefix, FormatParserStorage::getInstance().getTypeFromName(*iter)); } if (write_status) LOG(1, "All configs written."); else LOG(1, "Some config writing failed."); // store force index reference file { SortIndex_t SortIndex; BondFragments.StoreForcesFile(prefix, SortIndex); } // store keysets file TotalGraph.StoreKeySetFile(prefix); // store Hydrogen saturation correction file BondFragments.AddHydrogenCorrection(prefix); // restore orbital and Stop values //CalculateOrbitals(*configuration); } /** Internal helper to create from each keyset a molecule * */ void ExportGraph_ToFiles::prepareMolecule() { size_t count = 0; for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) { KeySet test = (*runner).first; LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "."); BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig())); ++count; } LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size() << " fragments generated from the keysets."); } /** Stores a fragment from \a KeySet into \a molecule. * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete * molecule and adds missing hydrogen where bonds were cut. * \param &Leaflet pointer to KeySet structure * \param IsAngstroem whether we have Ansgtroem or bohrradius * \return pointer to constructed molecule */ molecule * ExportGraph_ToFiles::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem) { Info info(__func__); ListOfLocalAtoms_t SonList; molecule *Leaf = World::getInstance().createMolecule(); StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList); // create the bonds between all: Make it an induced subgraph and add hydrogen // LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation)."); CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem); //Leaflet->Leaf->ScanForPeriodicCorrection(out); return Leaf; } /** Initializes some value for putting fragment of \a *mol into \a *Leaf. * \param *Leaf fragment molecule * \param &Leaflet pointer to KeySet structure * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol * \return number of atoms in fragment */ int ExportGraph_ToFiles::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList) { atom *FatherOfRunner = NULL; // first create the minimal set of atoms from the KeySet World &world = World::getInstance(); int size = 0; for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) { FatherOfRunner = world.getAtom(AtomById(*runner)); // find the id SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) ); size++; } return size; } /** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially). * \param *Leaf fragment molecule * \param IsAngstroem whether we have Ansgtroem or bohrradius * \param SonList list which atom of \a *Leaf is another atom's son */ void ExportGraph_ToFiles::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem) { bool LonelyFlag = false; atom *OtherFather = NULL; atom *FatherOfRunner = NULL; // we increment the iter just before skipping the hydrogen // as we use AddBond, we cannot have a const_iterator here for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) { LonelyFlag = true; FatherOfRunner = (*iter)->father; ASSERT(FatherOfRunner,"Atom without father found"); if (SonList.find(FatherOfRunner->getNr()) != SonList.end()) { // check if this, our father, is present in list // create all bonds const BondList& ListOfBonds = FatherOfRunner->getListOfBonds(); for (BondList::const_iterator BondRunner = ListOfBonds.begin(); BondRunner != ListOfBonds.end(); ++BondRunner) { OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner); if (SonList.find(OtherFather->getNr()) != SonList.end()) { // LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()] // << " is bound to " << *OtherFather << ", whose son is " // << *SonList[OtherFather->getNr()] << "."); if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba) std::stringstream output; // output << "ACCEPT: Adding Bond: " output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->BondDegree); // LOG(3, output.str()); //NumBonds[(*iter)->getNr()]++; } else { // LOG(3, "REJECY: Not adding bond, labels in wrong order."); } LonelyFlag = false; } else { // LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()] // << " is bound to " << *OtherFather << ", who has no son in this fragment molecule."); if (saturation == DoSaturate) { // LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between."); if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem)) exit(1); } //NumBonds[(*iter)->getNr()] += Binder->BondDegree; } } } else { ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!"); } if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) { LOG(0, **iter << "has got bonds only to hydrogens!"); } ++iter; if (saturation == DoSaturate) { while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen iter++; } } } }