[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 |
|
---|
[97ebf8] | 24 | /*
|
---|
| 25 | * FragmentationAction.cpp
|
---|
| 26 | *
|
---|
| 27 | * Created on: May 9, 2010
|
---|
| 28 | * Author: heber
|
---|
| 29 | */
|
---|
| 30 |
|
---|
[bf3817] | 31 | // include config.h
|
---|
| 32 | #ifdef HAVE_CONFIG_H
|
---|
| 33 | #include <config.h>
|
---|
| 34 | #endif
|
---|
| 35 |
|
---|
[ad011c] | 36 | #include "CodePatterns/MemDebug.hpp"
|
---|
[112b09] | 37 |
|
---|
[6f0841] | 38 | #include "Atom/atom.hpp"
|
---|
[c7c615] | 39 | #include "CodePatterns/IteratorAdaptors.hpp"
|
---|
[ad011c] | 40 | #include "CodePatterns/Log.hpp"
|
---|
[79d0b9] | 41 | #include "Descriptors/AtomSelectionDescriptor.hpp"
|
---|
[270bdf] | 42 | #include "Descriptors/MoleculeIdDescriptor.hpp"
|
---|
[d9dbef] | 43 | #include "Fragmentation/Exporters/ExportGraph_ToAtomFragments.hpp"
|
---|
[ca8bea] | 44 | #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
|
---|
[ac9ca4] | 45 | #include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
|
---|
[5d5550] | 46 | #include "Fragmentation/Exporters/SaturatedBond.hpp"
|
---|
[98a293b] | 47 | #include "Fragmentation/Exporters/SaturatedFragment.hpp"
|
---|
[5d5550] | 48 | #include "Fragmentation/Exporters/SaturationDistanceMaximizer.hpp"
|
---|
[246e13] | 49 | #include "Fragmentation/Fragmentation.hpp"
|
---|
[b4f72c] | 50 | #include "Fragmentation/Graph.hpp"
|
---|
[07a47e] | 51 | #include "Fragmentation/HydrogenSaturation_enum.hpp"
|
---|
[13c5c1] | 52 | #include "Fragmentation/Interfragmenter.hpp"
|
---|
[adb51ab] | 53 | #include "Fragmentation/KeySetsContainer.hpp"
|
---|
| 54 | #include "Fragmentation/Summation/Containers/FragmentationResultContainer.hpp"
|
---|
[0fad93] | 55 | #include "Graph/AdjacencyList.hpp"
|
---|
[79d0b9] | 56 | #include "Graph/BondGraph.hpp"
|
---|
[fe0cb8] | 57 | #include "Graph/CyclicStructureAnalysis.hpp"
|
---|
[49c059] | 58 | #include "Graph/DepthFirstSearchAnalysis.hpp"
|
---|
[9511c7] | 59 | #include "Helpers/defs.hpp"
|
---|
[97ebf8] | 60 | #include "molecule.hpp"
|
---|
| 61 | #include "World.hpp"
|
---|
| 62 |
|
---|
[c7c615] | 63 | #include <boost/shared_ptr.hpp>
|
---|
[3aa8a5] | 64 | #include <boost/filesystem.hpp>
|
---|
[c7c615] | 65 | #include <algorithm>
|
---|
[97ebf8] | 66 | #include <iostream>
|
---|
[2a0eb0] | 67 | #include <map>
|
---|
[97ebf8] | 68 | #include <string>
|
---|
[2a0eb0] | 69 | #include <vector>
|
---|
[97ebf8] | 70 |
|
---|
[1fd675] | 71 | #include "Actions/FragmentationAction/FragmentationAction.hpp"
|
---|
[70d9b9] | 72 |
|
---|
[ce7fdc] | 73 | using namespace MoleCuilder;
|
---|
| 74 |
|
---|
[1fd675] | 75 | // and construct the stuff
|
---|
| 76 | #include "FragmentationAction.def"
|
---|
| 77 | #include "Action_impl_pre.hpp"
|
---|
| 78 | /** =========== define the function ====================== */
|
---|
[b5b01e] | 79 | ActionState::ptr FragmentationFragmentationAction::performCall() {
|
---|
[e4b5de] | 80 | clock_t start,end;
|
---|
[2a0eb0] | 81 | int ExitFlag = -1;
|
---|
| 82 | World &world = World::getInstance();
|
---|
[e4b5de] | 83 |
|
---|
[2a0eb0] | 84 | // inform about used parameters
|
---|
[13c5c1] | 85 | LOG(0, "STATUS: Fragmenting molecular system with current connection matrix up to "
|
---|
[bae7bc] | 86 | << params.order.get() << " order. ");
|
---|
| 87 | if (params.types.get().size() != 0)
|
---|
| 88 | LOG(0, "STATUS: Fragment files begin with "
|
---|
| 89 | << params.prefix.get() << " and are stored as: "
|
---|
| 90 | << params.types.get() << "." << std::endl);
|
---|
[99b0dc] | 91 |
|
---|
[2a0eb0] | 92 | // check for selected atoms
|
---|
| 93 | if (world.beginAtomSelection() == world.endAtomSelection()) {
|
---|
[26b4d62] | 94 | STATUS("There are no atoms selected for fragmentation.");
|
---|
[2a0eb0] | 95 | return Action::failure;
|
---|
| 96 | }
|
---|
| 97 |
|
---|
| 98 | // go through all atoms, note down their molecules and group them
|
---|
[270bdf] | 99 | typedef std::multimap<const molecule *, atom *> clusters_t;
|
---|
[3aa8a5] | 100 | typedef std::vector<atomId_t> atomids_t;
|
---|
| 101 | atomids_t atomids;
|
---|
[2a0eb0] | 102 | clusters_t clusters;
|
---|
| 103 | for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
|
---|
| 104 | iter != world.endAtomSelection(); ++iter) {
|
---|
| 105 | clusters.insert( std::make_pair(iter->second->getMolecule(), iter->second) );
|
---|
[3aa8a5] | 106 | atomids.push_back(iter->second->getId());
|
---|
[2a0eb0] | 107 | }
|
---|
[c7c615] | 108 | {
|
---|
[270bdf] | 109 | std::vector<const molecule *> molecules;
|
---|
[c7c615] | 110 | molecules.insert( molecules.end(), MapKeyIterator<clusters_t::const_iterator>(clusters.begin()),
|
---|
| 111 | MapKeyIterator<clusters_t::const_iterator>(clusters.end()) );
|
---|
| 112 | molecules.erase( std::unique(molecules.begin(), molecules.end()), molecules.end() );
|
---|
| 113 | LOG(1, "INFO: There are " << molecules.size() << " molecules to consider.");
|
---|
| 114 | }
|
---|
[2a0eb0] | 115 |
|
---|
[9511c7] | 116 | // parse in Adjacency file
|
---|
[3aa8a5] | 117 | boost::shared_ptr<AdjacencyList> FileChecker;
|
---|
| 118 | boost::filesystem::path filename(params.prefix.get() + std::string(ADJACENCYFILE));
|
---|
| 119 | if (boost::filesystem::exists(filename) && boost::filesystem::is_regular_file(filename)) {
|
---|
| 120 | std::ifstream File;
|
---|
| 121 | File.open(filename.string().c_str(), ios::out);
|
---|
| 122 | FileChecker.reset(new AdjacencyList(File));
|
---|
| 123 | File.close();
|
---|
| 124 | } else {
|
---|
| 125 | LOG(1, "INFO: Could not open default adjacency file " << filename.string() << ".");
|
---|
| 126 | FileChecker.reset(new AdjacencyList);
|
---|
| 127 | }
|
---|
[9511c7] | 128 |
|
---|
[79d0b9] | 129 | // make sure bond degree is correct
|
---|
| 130 | {
|
---|
| 131 | BondGraph *BG = World::getInstance().getBondGraph();
|
---|
| 132 | World::AtomComposite Set = World::getInstance().getAllAtoms(AtomsBySelection());
|
---|
[0763ce] | 133 | // check whether bond graph is correct
|
---|
| 134 | if (!BG->checkBondDegree(Set))
|
---|
| 135 | BG->CorrectBondDegree(Set);
|
---|
| 136 | else
|
---|
| 137 | LOG(1, "INFO: Bond degrees all valid, not correcting.");
|
---|
[79d0b9] | 138 | }
|
---|
| 139 |
|
---|
[bfbd4a] | 140 | // we parse in the keysets from last time if present
|
---|
| 141 | Graph StoredGraph;
|
---|
| 142 | StoredGraph.ParseKeySetFile(params.prefix.get());
|
---|
| 143 |
|
---|
[2a0eb0] | 144 | start = clock();
|
---|
| 145 | // go through all keys (i.e. all molecules)
|
---|
| 146 | clusters_t::const_iterator advanceiter;
|
---|
[b4f72c] | 147 | Graph TotalGraph;
|
---|
| 148 | int keysetcounter = 0;
|
---|
[2a0eb0] | 149 | for (clusters_t::const_iterator iter = clusters.begin();
|
---|
| 150 | iter != clusters.end();
|
---|
| 151 | iter = advanceiter) {
|
---|
| 152 | // get iterator to past last atom in this molecule
|
---|
[270bdf] | 153 | const molecule * mol = iter->first;
|
---|
[2a0eb0] | 154 | advanceiter = clusters.upper_bound(mol);
|
---|
| 155 |
|
---|
| 156 | // copy molecule's atoms' ids as parameters to Fragmentation's AtomMask
|
---|
| 157 | std::vector<atomId_t> mols_atomids;
|
---|
| 158 | std::transform(iter, advanceiter, std::back_inserter(mols_atomids),
|
---|
| 159 | boost::bind( &atom::getNr,
|
---|
| 160 | boost::bind( &clusters_t::value_type::second, _1 ))
|
---|
| 161 | );
|
---|
| 162 | LOG(2, "INFO: Fragmenting in molecule " << mol->getName() << " in " << clusters.count(mol)
|
---|
| 163 | << " atoms, out of " << mol->getAtomCount() << ".");
|
---|
[9291d04] | 164 | const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
|
---|
[270bdf] | 165 | molecule * non_const_mol = World::getInstance().getMolecule(MoleculeById(mol->getId()));
|
---|
| 166 | Fragmentation Fragmenter(non_const_mol, *FileChecker, treatment);
|
---|
[2a0eb0] | 167 |
|
---|
| 168 | // perform fragmentation
|
---|
| 169 | LOG(0, std::endl << " ========== Fragmentation of molecule " << mol->getName() << " ========================= ");
|
---|
| 170 | {
|
---|
[bfbd4a] | 171 | Graph StoredLocalGraph(StoredGraph.getLocalGraph(mol));
|
---|
[fe0cb8] | 172 | const int tempFlag = Fragmenter.FragmentMolecule(mols_atomids, params.order.get(), params.prefix.get(), StoredLocalGraph);
|
---|
[2a0eb0] | 173 | if ((ExitFlag == 2) && (tempFlag != 2))
|
---|
| 174 | ExitFlag = tempFlag; // if there is one molecule that needs further fragmentation, it overrides others
|
---|
| 175 | if (ExitFlag == -1)
|
---|
| 176 | ExitFlag = tempFlag; // if we are the first, we set the standard
|
---|
[e4b5de] | 177 | }
|
---|
[569e42] | 178 | if (TotalGraph.empty()) {
|
---|
| 179 | TotalGraph = Fragmenter.getGraph();
|
---|
| 180 | keysetcounter = TotalGraph.size();
|
---|
| 181 | } else
|
---|
| 182 | TotalGraph.InsertGraph(Fragmenter.getGraph(), keysetcounter);
|
---|
[ca8bea] | 183 |
|
---|
[97ebf8] | 184 | }
|
---|
[fe0cb8] | 185 | // add full cycles if desired
|
---|
| 186 | if (params.DoCyclesFull.get()) {
|
---|
| 187 | // get the BackEdgeStack from somewhere
|
---|
| 188 | DepthFirstSearchAnalysis DFS;
|
---|
| 189 | DFS();
|
---|
| 190 | std::deque<bond::ptr> BackEdgeStack = DFS.getBackEdgeStack();
|
---|
| 191 | // then we analyse the cycles and get them
|
---|
| 192 | CyclicStructureAnalysis CycleAnalysis(params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen);
|
---|
| 193 | CycleAnalysis(&BackEdgeStack);
|
---|
| 194 | CyclicStructureAnalysis::cycles_t cycles = CycleAnalysis.getAllCycles();
|
---|
[adb51ab] | 195 | // sort them according to KeySet::operator<()
|
---|
| 196 | std::sort(cycles.begin(), cycles.end());
|
---|
| 197 | // store all found cycles to file
|
---|
| 198 | {
|
---|
| 199 | boost::filesystem::path filename(params.prefix.get() + std::string(CYCLEKEYSETFILE));
|
---|
| 200 | std::ofstream File;
|
---|
| 201 | LOG(1, "INFO: Storing cycle keysets to " << filename.string() << ".");
|
---|
| 202 | File.open(filename.string().c_str(), ios::out);
|
---|
| 203 | for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
|
---|
| 204 | iter != cycles.end(); ++iter) {
|
---|
| 205 | for (CyclicStructureAnalysis::cycle_t::const_iterator cycleiter = (*iter).begin();
|
---|
| 206 | cycleiter != (*iter).end(); ++cycleiter) {
|
---|
| 207 | File << *cycleiter << "\t";
|
---|
| 208 | }
|
---|
| 209 | File << "\n";
|
---|
| 210 | }
|
---|
| 211 | File.close();
|
---|
| 212 | }
|
---|
| 213 | // ... and to result container
|
---|
| 214 | {
|
---|
| 215 | KeySetsContainer cyclekeys;
|
---|
| 216 | for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
|
---|
| 217 | iter != cycles.end(); ++iter) {
|
---|
| 218 | const CyclicStructureAnalysis::cycle_t &cycle = *iter;
|
---|
| 219 | const size_t order = cycle.size();
|
---|
| 220 | KeySetsContainer::IntVector temp_cycle(cycle.begin(), cycle.end());
|
---|
| 221 | cyclekeys.insert(temp_cycle, order);
|
---|
| 222 | }
|
---|
| 223 | FragmentationResultContainer::getInstance().addCycles(cyclekeys);
|
---|
| 224 | }
|
---|
[fe0cb8] | 225 | // Create graph and insert into TotalGraph
|
---|
[adb51ab] | 226 | LOG(0, "STATUS: Adding " << cycles.size() << " cycles.");
|
---|
[fe0cb8] | 227 | {
|
---|
| 228 | Graph CycleGraph;
|
---|
| 229 | for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
|
---|
| 230 | iter != cycles.end(); ++iter) {
|
---|
| 231 | const CyclicStructureAnalysis::cycle_t ¤tcycle = *iter;
|
---|
| 232 | LOG(2, "INFO: Inserting cycle " << currentcycle << ".");
|
---|
| 233 | #ifndef NDEBUG
|
---|
| 234 | std::pair< Graph::iterator, bool > inserter =
|
---|
| 235 | #endif
|
---|
| 236 | CycleGraph.insert( std::make_pair(currentcycle, NumberValuePair(1,1.)) );
|
---|
| 237 | ASSERT( inserter.second,
|
---|
| 238 | "FragmentationFragmentationAction::performCall() - keyset "
|
---|
| 239 | +toString(currentcycle)+" inserted twice into CycleGraph.");
|
---|
| 240 | }
|
---|
| 241 | TotalGraph.InsertGraph(CycleGraph, keysetcounter);
|
---|
| 242 | }
|
---|
| 243 | }
|
---|
| 244 |
|
---|
[569e42] | 245 | LOG(0, "STATUS: There are " << TotalGraph.size() << " fragments.");
|
---|
[3aa8a5] | 246 |
|
---|
[321470] | 247 | {
|
---|
| 248 | // remove OrderAtSite file
|
---|
| 249 | std::string line;
|
---|
| 250 | std::ofstream file;
|
---|
| 251 | line = params.prefix.get() + ORDERATSITEFILE;
|
---|
| 252 | file.open(line.c_str(), std::ofstream::out | std::ofstream::trunc);
|
---|
| 253 | file << "";
|
---|
| 254 | file.close();
|
---|
| 255 | }
|
---|
| 256 |
|
---|
[d410e25] | 257 | // store graph internally
|
---|
| 258 | AtomFragmentsMap &atomfragments = AtomFragmentsMap::getInstance();
|
---|
| 259 | atomfragments.clear();
|
---|
| 260 | atomfragments.insert(TotalGraph);
|
---|
| 261 |
|
---|
[13c5c1] | 262 | // now add interfragments
|
---|
| 263 | if (params.InterOrder.get() != 0) {
|
---|
| 264 | LOG(0, "STATUS: Putting fragments together up to order "
|
---|
| 265 | << params.InterOrder.get() << " and distance of "
|
---|
| 266 | << params.distance.get() << ".");
|
---|
[cee9e8] | 267 | const enum HydrogenTreatment treatment =
|
---|
| 268 | params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
|
---|
| 269 | const double UpperBound = std::max(10., params.distance.get());
|
---|
[d410e25] | 270 | Interfragmenter fragmenter;
|
---|
[cee9e8] | 271 |
|
---|
| 272 | // check the largest Rcut that causes no additional inter-fragments
|
---|
| 273 | const double Min_Rcut =
|
---|
| 274 | fragmenter.findLargestCutoff(params.InterOrder.get(), UpperBound, treatment);
|
---|
| 275 |
|
---|
[92232f] | 276 | // if we smear out electronic charges, warn when non-overlapping criterion does not hold
|
---|
| 277 | if (params.InterOrder.get() < Min_Rcut)
|
---|
| 278 | ELOG(2, "Inter-order is too low to cause any additional fragments.");
|
---|
| 279 |
|
---|
[cee9e8] | 280 | // then add fragments
|
---|
[d410e25] | 281 | fragmenter(TotalGraph, params.InterOrder.get(), params.distance.get(), treatment);
|
---|
[cee9e8] | 282 |
|
---|
[13c5c1] | 283 | LOG(0, "STATUS: There are now " << TotalGraph.size() << " fragments after interfragmenting.");
|
---|
| 284 | }
|
---|
[d410e25] | 285 | // TODO: When insert only adds and replaces if already present, no clear needed
|
---|
[3004d2] | 286 | atomfragments.clear();
|
---|
| 287 | atomfragments.insert(TotalGraph);
|
---|
| 288 |
|
---|
[bfbd4a] | 289 | // store keysets to file
|
---|
| 290 | {
|
---|
| 291 | TotalGraph.StoreKeySetFile(params.prefix.get());
|
---|
| 292 | }
|
---|
| 293 |
|
---|
[98a293b] | 294 | // create global saturation positions map
|
---|
| 295 | SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
|
---|
[5d5550] | 296 | {
|
---|
| 297 | // go through each atom
|
---|
| 298 | for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
|
---|
| 299 | iter != world.endAtomSelection(); ++iter) {
|
---|
| 300 | const atom * const _atom = iter->second;
|
---|
| 301 |
|
---|
| 302 | // skip hydrogens if treated special
|
---|
| 303 | const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
|
---|
| 304 | if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {
|
---|
| 305 | LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);
|
---|
| 306 | continue;
|
---|
| 307 | }
|
---|
| 308 |
|
---|
| 309 | // get the valence
|
---|
| 310 | unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
|
---|
| 311 | LOG(3, "DEBUG: There are " << NumberOfPoints
|
---|
| 312 | << " places to fill in in total for this atom " << *_atom << ".");
|
---|
| 313 |
|
---|
| 314 | // check whether there are any bonds with degree larger than 1
|
---|
| 315 | unsigned int SumOfDegrees = 0;
|
---|
| 316 | bool PresentHigherBonds = false;
|
---|
| 317 | const BondList &bondlist = _atom->getListOfBonds();
|
---|
| 318 | for (BondList::const_iterator bonditer = bondlist.begin();
|
---|
| 319 | bonditer != bondlist.end(); ++bonditer) {
|
---|
| 320 | SumOfDegrees += (*bonditer)->getDegree();
|
---|
| 321 | PresentHigherBonds |= (*bonditer)->getDegree() > 1;
|
---|
| 322 | }
|
---|
| 323 |
|
---|
| 324 | // check whether there are alphas to maximize the hydrogens distances
|
---|
| 325 | SaturationDistanceMaximizer::position_bins_t position_bins;
|
---|
| 326 | {
|
---|
| 327 | // gather all bonds and convert to SaturatedBonds
|
---|
| 328 | SaturationDistanceMaximizer::PositionContainers_t CutBonds;
|
---|
| 329 | for (BondList::const_iterator bonditer = bondlist.begin();
|
---|
| 330 | bonditer != bondlist.end(); ++bonditer) {
|
---|
| 331 | CutBonds.push_back(
|
---|
| 332 | SaturatedBond::ptr(new SaturatedBond(*(bonditer->get()), *_atom) )
|
---|
| 333 | );
|
---|
| 334 | }
|
---|
| 335 | SaturationDistanceMaximizer maximizer(CutBonds);
|
---|
| 336 | if (PresentHigherBonds) {
|
---|
| 337 | // then find best alphas
|
---|
| 338 | maximizer();
|
---|
| 339 | } else {
|
---|
| 340 | // if no higher order bonds, we simply gather the scaled positions
|
---|
| 341 | }
|
---|
| 342 | position_bins = maximizer.getAllPositionBins();
|
---|
| 343 | LOG(4, "DEBUG: Positions for atom " << *_atom << " are " << position_bins);
|
---|
| 344 | }
|
---|
| 345 |
|
---|
| 346 | // convert into the desired entry in the map
|
---|
| 347 | SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;
|
---|
| 348 | {
|
---|
| 349 | BondList::const_iterator bonditer = bondlist.begin();
|
---|
| 350 | SaturationDistanceMaximizer::position_bins_t::const_iterator biniter =
|
---|
| 351 | position_bins.begin();
|
---|
| 352 |
|
---|
| 353 | for (;bonditer != bondlist.end(); ++bonditer, ++biniter) {
|
---|
| 354 | const atom * const OtherAtom = (*bonditer)->GetOtherAtom(_atom);
|
---|
| 355 | std::pair<
|
---|
| 356 | SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,
|
---|
| 357 | bool
|
---|
| 358 | > inserter;
|
---|
| 359 | // check whether we treat hydrogen special
|
---|
| 360 | if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {
|
---|
| 361 | // if hydrogen, forget rescaled position and use original one
|
---|
| 362 | inserter =
|
---|
| 363 | positions_per_neighbor.insert(
|
---|
| 364 | std::make_pair(
|
---|
| 365 | OtherAtom->getId(),
|
---|
| 366 | SaturatedFragment::SaturationsPositions_t(
|
---|
| 367 | 1, OtherAtom->getPosition() - _atom->getPosition())
|
---|
| 368 | )
|
---|
| 369 | );
|
---|
| 370 | } else {
|
---|
| 371 | inserter =
|
---|
| 372 | positions_per_neighbor.insert(
|
---|
| 373 | std::make_pair(
|
---|
| 374 | OtherAtom->getId(),
|
---|
| 375 | SaturatedFragment::SaturationsPositions_t(
|
---|
| 376 | biniter->begin(),
|
---|
| 377 | biniter->end())
|
---|
| 378 | )
|
---|
| 379 | );
|
---|
| 380 | }
|
---|
| 381 | // if already pressent, add to this present list
|
---|
| 382 | ASSERT (inserter.second,
|
---|
| 383 | "FragmentationAction::performCall() - other atom "
|
---|
| 384 | +toString(*OtherAtom)+" already present?");
|
---|
| 385 | }
|
---|
| 386 | // bonditer follows nicely
|
---|
| 387 | ASSERT( biniter == position_bins.end(),
|
---|
| 388 | "FragmentationAction::performCall() - biniter is out of step, it still points at bond "
|
---|
| 389 | +toString(*biniter)+".");
|
---|
| 390 | }
|
---|
| 391 | // and insert
|
---|
| 392 | globalsaturationpositions.insert(
|
---|
| 393 | std::make_pair( _atom->getId(),
|
---|
| 394 | positions_per_neighbor
|
---|
| 395 | ));
|
---|
| 396 | }
|
---|
| 397 | }
|
---|
[98a293b] | 398 |
|
---|
[3aa8a5] | 399 | {
|
---|
[9291d04] | 400 | const enum HydrogenSaturation saturation = params.DoSaturation.get() ? DoSaturate : DontSaturate;
|
---|
[276ac6] | 401 | const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
|
---|
[ac9ca4] | 402 | if (params.types.get().size() != 0) {
|
---|
| 403 | // store molecule's fragment to file
|
---|
[98a293b] | 404 | ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
|
---|
[ac9ca4] | 405 | exporter.setPrefix(params.prefix.get());
|
---|
| 406 | exporter.setOutputTypes(params.types.get());
|
---|
| 407 | exporter();
|
---|
| 408 | } else {
|
---|
| 409 | // store molecule's fragment in FragmentJobQueue
|
---|
[98a293b] | 410 | ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
|
---|
[ac9ca4] | 411 | exporter.setLevel(params.level.get());
|
---|
| 412 | exporter();
|
---|
| 413 | }
|
---|
[d9dbef] | 414 | ExportGraph_ToAtomFragments exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
|
---|
| 415 | exporter();
|
---|
| 416 | }
|
---|
| 417 | if (!AtomFragmentsMap::getInstance().checkCompleteness()) {
|
---|
| 418 | ELOG(0, "Something went wrong with placing keysets in AtomFragmentsMap.");
|
---|
| 419 | return Action::failure;
|
---|
[3aa8a5] | 420 | }
|
---|
| 421 |
|
---|
| 422 | // store Adjacency to file
|
---|
| 423 | {
|
---|
| 424 | std::string filename = params.prefix.get() + ADJACENCYFILE;
|
---|
| 425 | std::ofstream AdjacencyFile;
|
---|
| 426 | AdjacencyFile.open(filename.c_str(), ios::out);
|
---|
| 427 | AdjacencyList adjacency(atomids);
|
---|
| 428 | adjacency.StoreToFile(AdjacencyFile);
|
---|
| 429 | AdjacencyFile.close();
|
---|
| 430 | }
|
---|
| 431 |
|
---|
[2a0eb0] | 432 | World::getInstance().setExitFlag(ExitFlag);
|
---|
| 433 | end = clock();
|
---|
| 434 | LOG(0, "STATUS: Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s.");
|
---|
| 435 |
|
---|
[70d9b9] | 436 | return Action::success;
|
---|
[97ebf8] | 437 | }
|
---|
| 438 |
|
---|
[b5b01e] | 439 | ActionState::ptr FragmentationFragmentationAction::performUndo(ActionState::ptr _state) {
|
---|
[70d9b9] | 440 | return Action::success;
|
---|
[97ebf8] | 441 | }
|
---|
| 442 |
|
---|
[b5b01e] | 443 | ActionState::ptr FragmentationFragmentationAction::performRedo(ActionState::ptr _state){
|
---|
[70d9b9] | 444 | return Action::success;
|
---|
[97ebf8] | 445 | }
|
---|
| 446 |
|
---|
| 447 | bool FragmentationFragmentationAction::canUndo() {
|
---|
[70d9b9] | 448 | return true;
|
---|
[97ebf8] | 449 | }
|
---|
| 450 |
|
---|
| 451 | bool FragmentationFragmentationAction::shouldUndo() {
|
---|
[70d9b9] | 452 | return true;
|
---|
[97ebf8] | 453 | }
|
---|
[1fd675] | 454 | /** =========== end of function ====================== */
|
---|