source: src/Actions/FragmentationAction/FragmentationAction.cpp@ 56c55f3

FragmentMolecule_checks_bonddegrees Gui_Fixes
Last change on this file since 56c55f3 was 56c55f3, checked in by Frederik Heber <heber@…>, 9 years ago

Added "parse-state-files" option to fragment-molecule which is off by default.

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