source: src/Actions/FragmentationAction/FragmentationAction.cpp@ ce70d25

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes stable
Last change on this file since ce70d25 was 0f5956, checked in by Frederik Heber <heber@…>, 8 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.8 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 }
142
143 // we parse in the keysets from last time if present
144 Graph StoredGraph;
145 if (params.ParseStateFiles.get()) {
146 StoredGraph.ParseKeySetFile(params.prefix.get());
147 // check parsed graph against the set of atoms
148 {
149 AdaptivityMap *amap = StoredGraph.GraphToAdaptivityMap();
150 bool status = true;
151 for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
152 iter != world.endAtomSelection(); ++iter) {
153 const atomId_t atomid = iter->second->getId();
154 // skip hydrogens in check if saturation is turned on
155 if ((iter->second->getType()->getAtomicNumber() != 1)
156 || (!params.DoSaturation.get())) {
157 if (amap->count(atomid) == 0) {
158 ELOG(1, "Atom #" << atomid << " not contained in KeySet file. ");
159 status = false;
160 break;
161 }
162 } else if (amap->count(atomid) != 0) {
163 ELOG(1, "Atom #" << atomid << " in KeySet file is a hydrogen, but is now excluded. ");
164 status = false;
165 break;
166 }
167 }
168 delete amap;
169
170 if (!status) {
171 ELOG(1, "KeySetsFile seems to contain leftover from an old fragmentation, hence not using file.");
172 StoredGraph.clear();
173 }
174 }
175 }
176
177 start = clock();
178 // go through all keys (i.e. all molecules)
179 clusters_t::const_iterator advanceiter;
180 Graph TotalGraph;
181 int keysetcounter = 0;
182 for (clusters_t::const_iterator iter = clusters.begin();
183 iter != clusters.end();
184 iter = advanceiter) {
185 // get iterator to past last atom in this molecule
186 const molecule * mol = iter->first;
187 advanceiter = clusters.upper_bound(mol);
188
189 // copy molecule's atoms' ids as parameters to Fragmentation's AtomMask
190 std::vector<atomId_t> mols_atomids;
191 std::transform(iter, advanceiter, std::back_inserter(mols_atomids),
192 boost::bind( &atom::getNr,
193 boost::bind( &clusters_t::value_type::second, _1 ))
194 );
195 LOG(2, "INFO: Fragmenting in molecule " << mol->getName() << " in " << clusters.count(mol)
196 << " atoms, out of " << mol->getAtomCount() << ".");
197 const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
198 molecule * non_const_mol = World::getInstance().getMolecule(MoleculeById(mol->getId()));
199 Fragmentation Fragmenter(non_const_mol, *FileChecker, treatment);
200
201 // perform fragmentation
202 LOG(0, std::endl << " ========== Fragmentation of molecule " << mol->getName() << " ========================= ");
203 {
204 Graph StoredLocalGraph(StoredGraph.getLocalGraph(mol));
205 const int tempFlag = Fragmenter.FragmentMolecule(mols_atomids, params.order.get(), params.prefix.get(), StoredLocalGraph, params.ParseStateFiles.get());
206 if ((ExitFlag == 2) && (tempFlag != 2))
207 ExitFlag = tempFlag; // if there is one molecule that needs further fragmentation, it overrides others
208 if (ExitFlag == -1)
209 ExitFlag = tempFlag; // if we are the first, we set the standard
210 }
211 if (TotalGraph.empty()) {
212 TotalGraph = Fragmenter.getGraph();
213 keysetcounter = TotalGraph.size();
214 } else
215 TotalGraph.InsertGraph(Fragmenter.getGraph(), keysetcounter);
216
217 }
218 // add full cycles if desired
219 if (params.DoCyclesFull.get()) {
220 // get the BackEdgeStack from somewhere
221 DepthFirstSearchAnalysis DFS;
222 DFS();
223 std::deque<bond::ptr> BackEdgeStack = DFS.getBackEdgeStack();
224 // then we analyse the cycles and get them
225 CyclicStructureAnalysis CycleAnalysis(params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen);
226 CycleAnalysis(&BackEdgeStack);
227 CyclicStructureAnalysis::cycles_t cycles = CycleAnalysis.getAllCycles();
228 // sort them according to KeySet::operator<()
229 std::sort(cycles.begin(), cycles.end());
230 // store all found cycles to file
231 {
232 boost::filesystem::path filename(params.prefix.get() + std::string(CYCLEKEYSETFILE));
233 std::ofstream File;
234 LOG(1, "INFO: Storing cycle keysets to " << filename.string() << ".");
235 File.open(filename.string().c_str(), ios::out);
236 for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
237 iter != cycles.end(); ++iter) {
238 for (CyclicStructureAnalysis::cycle_t::const_iterator cycleiter = (*iter).begin();
239 cycleiter != (*iter).end(); ++cycleiter) {
240 File << *cycleiter << "\t";
241 }
242 File << "\n";
243 }
244 File.close();
245 }
246 // ... and to result container
247 {
248 KeySetsContainer cyclekeys;
249 for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
250 iter != cycles.end(); ++iter) {
251 const CyclicStructureAnalysis::cycle_t &cycle = *iter;
252 const size_t order = cycle.size();
253 KeySetsContainer::IntVector temp_cycle(cycle.begin(), cycle.end());
254 cyclekeys.insert(temp_cycle, order);
255 }
256 FragmentationResultContainer::getInstance().addCycles(cyclekeys);
257 }
258 // Create graph and insert into TotalGraph
259 LOG(0, "STATUS: Adding " << cycles.size() << " cycles.");
260 {
261 Graph CycleGraph;
262 for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
263 iter != cycles.end(); ++iter) {
264 const CyclicStructureAnalysis::cycle_t &currentcycle = *iter;
265 LOG(2, "INFO: Inserting cycle " << currentcycle << ".");
266#ifndef NDEBUG
267 std::pair< Graph::iterator, bool > inserter =
268#endif
269 CycleGraph.insert( std::make_pair(currentcycle, NumberValuePair(1,1.)) );
270 ASSERT( inserter.second,
271 "FragmentationFragmentationAction::performCall() - keyset "
272 +toString(currentcycle)+" inserted twice into CycleGraph.");
273 }
274 TotalGraph.InsertGraph(CycleGraph, keysetcounter);
275 }
276 }
277
278 LOG(0, "STATUS: There are " << TotalGraph.size() << " fragments.");
279
280 {
281 // remove OrderAtSite file
282 std::string line;
283 std::ofstream file;
284 line = params.prefix.get() + ORDERATSITEFILE;
285 file.open(line.c_str(), std::ofstream::out | std::ofstream::trunc);
286 file << "";
287 file.close();
288 }
289
290 // store graph internally
291 AtomFragmentsMap &atomfragments = AtomFragmentsMap::getInstance();
292 atomfragments.clear();
293 atomfragments.insert(TotalGraph);
294
295 // now add interfragments
296 if (params.InterOrder.get() != 0) {
297 LOG(0, "STATUS: Putting fragments together up to order "
298 << params.InterOrder.get() << " and distance of "
299 << params.distance.get() << ".");
300 const enum HydrogenTreatment treatment =
301 params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
302 const double UpperBound = std::max(10., params.distance.get());
303 Interfragmenter fragmenter;
304
305 // check the largest Rcut that causes no additional inter-fragments
306 const double Min_Rcut =
307 fragmenter.findLargestCutoff(params.InterOrder.get(), UpperBound, treatment);
308
309 // if we smear out electronic charges, warn when non-overlapping criterion does not hold
310 if (params.InterOrder.get() < Min_Rcut)
311 ELOG(2, "Inter-order is too low to cause any additional fragments.");
312
313 // then add fragments
314 fragmenter(TotalGraph, params.InterOrder.get(), params.distance.get(), treatment);
315
316 LOG(0, "STATUS: There are now " << TotalGraph.size() << " fragments after interfragmenting.");
317 }
318 // TODO: When insert only adds and replaces if already present, no clear needed
319 atomfragments.clear();
320 atomfragments.insert(TotalGraph);
321
322 // store keysets to file
323 {
324 TotalGraph.StoreKeySetFile(params.prefix.get());
325 }
326
327 // create global saturation positions map
328 SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
329 {
330 // go through each atom
331 for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
332 iter != world.endAtomSelection(); ++iter) {
333 const atom * const _atom = iter->second;
334
335 // skip hydrogens if treated special
336 const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
337 if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {
338 LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);
339 continue;
340 }
341
342 // get the valence
343 unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
344 LOG(3, "DEBUG: There are " << NumberOfPoints
345 << " places to fill in in total for this atom " << *_atom << ".");
346
347 // check whether there are any bonds with degree larger than 1
348 unsigned int SumOfDegrees = 0;
349 bool PresentHigherBonds = false;
350 const BondList &bondlist = _atom->getListOfBonds();
351 for (BondList::const_iterator bonditer = bondlist.begin();
352 bonditer != bondlist.end(); ++bonditer) {
353 SumOfDegrees += (*bonditer)->getDegree();
354 PresentHigherBonds |= (*bonditer)->getDegree() > 1;
355 }
356
357 // check whether there are alphas to maximize the hydrogens distances
358 SaturationDistanceMaximizer::position_bins_t position_bins;
359 {
360 // gather all bonds and convert to SaturatedBonds
361 SaturationDistanceMaximizer::PositionContainers_t CutBonds;
362 for (BondList::const_iterator bonditer = bondlist.begin();
363 bonditer != bondlist.end(); ++bonditer) {
364 CutBonds.push_back(
365 SaturatedBond::ptr(new SaturatedBond(*(bonditer->get()), *_atom) )
366 );
367 }
368 SaturationDistanceMaximizer maximizer(CutBonds);
369 if (PresentHigherBonds) {
370 // then find best alphas
371 maximizer();
372 } else {
373 // if no higher order bonds, we simply gather the scaled positions
374 }
375 position_bins = maximizer.getAllPositionBins();
376 LOG(4, "DEBUG: Positions for atom " << *_atom << " are " << position_bins);
377 }
378
379 // convert into the desired entry in the map
380 SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;
381 {
382 BondList::const_iterator bonditer = bondlist.begin();
383 SaturationDistanceMaximizer::position_bins_t::const_iterator biniter =
384 position_bins.begin();
385
386 for (;bonditer != bondlist.end(); ++bonditer, ++biniter) {
387 const atom * const OtherAtom = (*bonditer)->GetOtherAtom(_atom);
388 std::pair<
389 SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,
390 bool
391 > inserter;
392 // check whether we treat hydrogen special
393 if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {
394 // if hydrogen, forget rescaled position and use original one
395 inserter =
396 positions_per_neighbor.insert(
397 std::make_pair(
398 OtherAtom->getId(),
399 SaturatedFragment::SaturationsPositions_t(
400 1, OtherAtom->getPosition() - _atom->getPosition())
401 )
402 );
403 } else {
404 inserter =
405 positions_per_neighbor.insert(
406 std::make_pair(
407 OtherAtom->getId(),
408 SaturatedFragment::SaturationsPositions_t(
409 biniter->begin(),
410 biniter->end())
411 )
412 );
413 }
414 // if already pressent, add to this present list
415 ASSERT (inserter.second,
416 "FragmentationAction::performCall() - other atom "
417 +toString(*OtherAtom)+" already present?");
418 }
419 // bonditer follows nicely
420 ASSERT( biniter == position_bins.end(),
421 "FragmentationAction::performCall() - biniter is out of step, it still points at bond "
422 +toString(*biniter)+".");
423 }
424 // and insert
425 globalsaturationpositions.insert(
426 std::make_pair( _atom->getId(),
427 positions_per_neighbor
428 ));
429 }
430 }
431
432 {
433 const enum HydrogenSaturation saturation = params.DoSaturation.get() ? DoSaturate : DontSaturate;
434 const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
435 if (params.types.get().size() != 0) {
436 // store molecule's fragment to file
437 ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
438 exporter.setPrefix(params.prefix.get());
439 exporter.setOutputTypes(params.types.get());
440 exporter();
441 } else {
442 // store molecule's fragment in FragmentJobQueue
443 ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
444 exporter.setLevel(params.level.get());
445 exporter();
446 }
447 ExportGraph_ToAtomFragments exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
448 exporter();
449 }
450 if (!AtomFragmentsMap::getInstance().checkCompleteness()) {
451 ELOG(0, "Something went wrong with placing keysets in AtomFragmentsMap.");
452 return Action::failure;
453 }
454
455 // store Adjacency to file
456 {
457 std::string filename = params.prefix.get() + ADJACENCYFILE;
458 std::ofstream AdjacencyFile;
459 AdjacencyFile.open(filename.c_str(), ios::out);
460 AdjacencyList adjacency(atomids);
461 adjacency.StoreToFile(AdjacencyFile);
462 AdjacencyFile.close();
463 }
464
465 World::getInstance().setExitFlag(ExitFlag);
466 end = clock();
467 LOG(0, "STATUS: Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s.");
468
469 return Action::success;
470}
471
472ActionState::ptr FragmentationFragmentationAction::performUndo(ActionState::ptr _state) {
473 return Action::success;
474}
475
476ActionState::ptr FragmentationFragmentationAction::performRedo(ActionState::ptr _state){
477 return Action::success;
478}
479
480bool FragmentationFragmentationAction::canUndo() {
481 return true;
482}
483
484bool FragmentationFragmentationAction::shouldUndo() {
485 return true;
486}
487/** =========== end of function ====================== */
Note: See TracBrowser for help on using the repository browser.