source: src/Actions/FragmentationAction/FragmentationAction.cpp@ 4a33f6

Last change on this file since 4a33f6 was 4a33f6, checked in by Frederik Heber <heber@…>, 10 years ago

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