/* * 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 . */ /* * FragmentationAction.cpp * * Created on: May 9, 2010 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "Atom/atom.hpp" #include "CodePatterns/IteratorAdaptors.hpp" #include "CodePatterns/Log.hpp" #include "Descriptors/AtomSelectionDescriptor.hpp" #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp" #ifdef HAVE_JOBMARKET #include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp" #endif #include "Fragmentation/Fragmentation.hpp" #include "Fragmentation/Graph.hpp" #include "Fragmentation/HydrogenSaturation_enum.hpp" #include "Graph/AdjacencyList.hpp" #include "Graph/BondGraph.hpp" #include "Graph/DepthFirstSearchAnalysis.hpp" #include "Helpers/defs.hpp" #include "molecule.hpp" #include "World.hpp" #include #include #include #include #include #include #include #include "Actions/FragmentationAction/FragmentationAction.hpp" using namespace MoleCuilder; // and construct the stuff #include "FragmentationAction.def" #include "Action_impl_pre.hpp" /** =========== define the function ====================== */ Action::state_ptr FragmentationFragmentationAction::performCall() { clock_t start,end; int ExitFlag = -1; World &world = World::getInstance(); // inform about used parameters LOG(0, "STATUS: Fragmenting molecular system with current connection matrix maximum bond distance " << params.distance.get() << " up to " << params.order.get() << " order. "); if (params.types.get().size() != 0) LOG(0, "STATUS: Fragment files begin with " << params.prefix.get() << " and are stored as: " << params.types.get() << "." << std::endl); // check for selected atoms if (world.beginAtomSelection() == world.endAtomSelection()) { ELOG(1, "There are not atoms selected for fragmentation."); return Action::failure; } // go through all atoms, note down their molecules and group them typedef std::multimap clusters_t; typedef std::vector atomids_t; atomids_t atomids; clusters_t clusters; for (World::AtomSelectionConstIterator iter = world.beginAtomSelection(); iter != world.endAtomSelection(); ++iter) { clusters.insert( std::make_pair(iter->second->getMolecule(), iter->second) ); atomids.push_back(iter->second->getId()); } { std::vector molecules; molecules.insert( molecules.end(), MapKeyIterator(clusters.begin()), MapKeyIterator(clusters.end()) ); molecules.erase( std::unique(molecules.begin(), molecules.end()), molecules.end() ); LOG(1, "INFO: There are " << molecules.size() << " molecules to consider."); } // parse in Adjacency file boost::shared_ptr FileChecker; boost::filesystem::path filename(params.prefix.get() + std::string(ADJACENCYFILE)); if (boost::filesystem::exists(filename) && boost::filesystem::is_regular_file(filename)) { std::ifstream File; File.open(filename.string().c_str(), ios::out); FileChecker.reset(new AdjacencyList(File)); File.close(); } else { LOG(1, "INFO: Could not open default adjacency file " << filename.string() << "."); FileChecker.reset(new AdjacencyList); } // make sure bond degree is correct { BondGraph *BG = World::getInstance().getBondGraph(); World::AtomComposite Set = World::getInstance().getAllAtoms(AtomsBySelection()); BG->CorrectBondDegree(Set); } // we require the current bond graph DepthFirstSearchAnalysis DFS; // we parse in the keysets from last time if present Graph StoredGraph; StoredGraph.ParseKeySetFile(params.prefix.get()); start = clock(); // go through all keys (i.e. all molecules) clusters_t::const_iterator advanceiter; Graph TotalGraph; int keysetcounter = 0; for (clusters_t::const_iterator iter = clusters.begin(); iter != clusters.end(); iter = advanceiter) { // get iterator to past last atom in this molecule molecule * mol = iter->first; advanceiter = clusters.upper_bound(mol); // copy molecule's atoms' ids as parameters to Fragmentation's AtomMask std::vector mols_atomids; std::transform(iter, advanceiter, std::back_inserter(mols_atomids), boost::bind( &atom::getNr, boost::bind( &clusters_t::value_type::second, _1 )) ); LOG(2, "INFO: Fragmenting in molecule " << mol->getName() << " in " << clusters.count(mol) << " atoms, out of " << mol->getAtomCount() << "."); const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen; Fragmentation Fragmenter(mol, *FileChecker, treatment); // perform fragmentation LOG(0, std::endl << " ========== Fragmentation of molecule " << mol->getName() << " ========================= "); { Graph StoredLocalGraph(StoredGraph.getLocalGraph(mol)); const int tempFlag = Fragmenter.FragmentMolecule(mols_atomids, params.order.get(), params.prefix.get(), DFS, StoredLocalGraph); if ((ExitFlag == 2) && (tempFlag != 2)) ExitFlag = tempFlag; // if there is one molecule that needs further fragmentation, it overrides others if (ExitFlag == -1) ExitFlag = tempFlag; // if we are the first, we set the standard } if (TotalGraph.empty()) { TotalGraph = Fragmenter.getGraph(); keysetcounter = TotalGraph.size(); } else TotalGraph.InsertGraph(Fragmenter.getGraph(), keysetcounter); } LOG(0, "STATUS: There are " << TotalGraph.size() << " fragments."); // store keysets to file { TotalGraph.StoreKeySetFile(params.prefix.get()); } { const enum HydrogenSaturation saturation = params.DoSaturation.get() ? DoSaturate : DontSaturate; const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen; if (params.types.get().size() != 0) { // store molecule's fragment to file ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation); exporter.setPrefix(params.prefix.get()); exporter.setOutputTypes(params.types.get()); exporter(); } else { #ifdef HAVE_JOBMARKET // store molecule's fragment in FragmentJobQueue ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation); exporter.setLevel(params.level.get()); exporter(); #else ELOG(1, "No output file types specified and JobMarket support is not compiled in."); return Action::failure; #endif } } // store Adjacency to file { std::string filename = params.prefix.get() + ADJACENCYFILE; std::ofstream AdjacencyFile; AdjacencyFile.open(filename.c_str(), ios::out); AdjacencyList adjacency(atomids); adjacency.StoreToFile(AdjacencyFile); AdjacencyFile.close(); } World::getInstance().setExitFlag(ExitFlag); end = clock(); LOG(0, "STATUS: Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s."); return Action::success; } Action::state_ptr FragmentationFragmentationAction::performUndo(Action::state_ptr _state) { return Action::success; } Action::state_ptr FragmentationFragmentationAction::performRedo(Action::state_ptr _state){ return Action::success; } bool FragmentationFragmentationAction::canUndo() { return true; } bool FragmentationFragmentationAction::shouldUndo() { return true; } /** =========== end of function ====================== */