source: src/Actions/FragmentationAction/FragmentationAction.cpp@ 70aeed

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing 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_PopActions 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 FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes 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 PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 70aeed was d410e25, checked in by Frederik Heber <heber@…>, 9 years ago

Interfragmenter now access singleton AtomFragmentsMap.

  • i.e. we do not instantiate the map ourselves but assume that it has been filled elsewhere. We only receive a Graph ("TotalGraph") in order to add the additional fragments to it.
  • FragmentationAction fills the AtomFragmentsMap already before calling interfragmenter.
  • Property mode set to 100644
File size: 17.2 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 FragmentationResultContainer::getInstance().addCycles(cyclekeys);
223 }
224 // Create graph and insert into TotalGraph
225 LOG(0, "STATUS: Adding " << cycles.size() << " cycles.");
226 {
227 Graph CycleGraph;
228 for (CyclicStructureAnalysis::cycles_t::const_iterator iter = cycles.begin();
229 iter != cycles.end(); ++iter) {
230 const CyclicStructureAnalysis::cycle_t &currentcycle = *iter;
231 LOG(2, "INFO: Inserting cycle " << currentcycle << ".");
232#ifndef NDEBUG
233 std::pair< Graph::iterator, bool > inserter =
234#endif
235 CycleGraph.insert( std::make_pair(currentcycle, NumberValuePair(1,1.)) );
236 ASSERT( inserter.second,
237 "FragmentationFragmentationAction::performCall() - keyset "
238 +toString(currentcycle)+" inserted twice into CycleGraph.");
239 }
240 TotalGraph.InsertGraph(CycleGraph, keysetcounter);
241 }
242 }
243
244 LOG(0, "STATUS: There are " << TotalGraph.size() << " fragments.");
245
246 {
247 // remove OrderAtSite file
248 std::string line;
249 std::ofstream file;
250 line = params.prefix.get() + ORDERATSITEFILE;
251 file.open(line.c_str(), std::ofstream::out | std::ofstream::trunc);
252 file << "";
253 file.close();
254 }
255
256 // store graph internally
257 AtomFragmentsMap &atomfragments = AtomFragmentsMap::getInstance();
258 atomfragments.clear();
259 atomfragments.insert(TotalGraph);
260
261 // now add interfragments
262 if (params.InterOrder.get() != 0) {
263 LOG(0, "STATUS: Putting fragments together up to order "
264 << params.InterOrder.get() << " and distance of "
265 << params.distance.get() << ".");
266 const enum HydrogenTreatment treatment =
267 params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
268 const double UpperBound = std::max(10., params.distance.get());
269 Interfragmenter fragmenter;
270
271 // check the largest Rcut that causes no additional inter-fragments
272 const double Min_Rcut =
273 fragmenter.findLargestCutoff(params.InterOrder.get(), UpperBound, treatment);
274
275 // if we smear out electronic charges, warn when non-overlapping criterion does not hold
276 if (params.InterOrder.get() < Min_Rcut)
277 ELOG(2, "Inter-order is too low to cause any additional fragments.");
278
279 // then add fragments
280 fragmenter(TotalGraph, params.InterOrder.get(), params.distance.get(), treatment);
281
282 LOG(0, "STATUS: There are now " << TotalGraph.size() << " fragments after interfragmenting.");
283 }
284 // TODO: When insert only adds and replaces if already present, no clear needed
285 atomfragments.clear();
286 atomfragments.insert(TotalGraph);
287
288 // store keysets to file
289 {
290 TotalGraph.StoreKeySetFile(params.prefix.get());
291 }
292
293 // create global saturation positions map
294 SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
295 {
296 // go through each atom
297 for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
298 iter != world.endAtomSelection(); ++iter) {
299 const atom * const _atom = iter->second;
300
301 // skip hydrogens if treated special
302 const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
303 if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {
304 LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);
305 continue;
306 }
307
308 // get the valence
309 unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
310 LOG(3, "DEBUG: There are " << NumberOfPoints
311 << " places to fill in in total for this atom " << *_atom << ".");
312
313 // check whether there are any bonds with degree larger than 1
314 unsigned int SumOfDegrees = 0;
315 bool PresentHigherBonds = false;
316 const BondList &bondlist = _atom->getListOfBonds();
317 for (BondList::const_iterator bonditer = bondlist.begin();
318 bonditer != bondlist.end(); ++bonditer) {
319 SumOfDegrees += (*bonditer)->getDegree();
320 PresentHigherBonds |= (*bonditer)->getDegree() > 1;
321 }
322
323 // check whether there are alphas to maximize the hydrogens distances
324 SaturationDistanceMaximizer::position_bins_t position_bins;
325 {
326 // gather all bonds and convert to SaturatedBonds
327 SaturationDistanceMaximizer::PositionContainers_t CutBonds;
328 for (BondList::const_iterator bonditer = bondlist.begin();
329 bonditer != bondlist.end(); ++bonditer) {
330 CutBonds.push_back(
331 SaturatedBond::ptr(new SaturatedBond(*(bonditer->get()), *_atom) )
332 );
333 }
334 SaturationDistanceMaximizer maximizer(CutBonds);
335 if (PresentHigherBonds) {
336 // then find best alphas
337 maximizer();
338 } else {
339 // if no higher order bonds, we simply gather the scaled positions
340 }
341 position_bins = maximizer.getAllPositionBins();
342 LOG(4, "DEBUG: Positions for atom " << *_atom << " are " << position_bins);
343 }
344
345 // convert into the desired entry in the map
346 SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;
347 {
348 BondList::const_iterator bonditer = bondlist.begin();
349 SaturationDistanceMaximizer::position_bins_t::const_iterator biniter =
350 position_bins.begin();
351
352 for (;bonditer != bondlist.end(); ++bonditer, ++biniter) {
353 const atom * const OtherAtom = (*bonditer)->GetOtherAtom(_atom);
354 std::pair<
355 SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,
356 bool
357 > inserter;
358 // check whether we treat hydrogen special
359 if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {
360 // if hydrogen, forget rescaled position and use original one
361 inserter =
362 positions_per_neighbor.insert(
363 std::make_pair(
364 OtherAtom->getId(),
365 SaturatedFragment::SaturationsPositions_t(
366 1, OtherAtom->getPosition() - _atom->getPosition())
367 )
368 );
369 } else {
370 inserter =
371 positions_per_neighbor.insert(
372 std::make_pair(
373 OtherAtom->getId(),
374 SaturatedFragment::SaturationsPositions_t(
375 biniter->begin(),
376 biniter->end())
377 )
378 );
379 }
380 // if already pressent, add to this present list
381 ASSERT (inserter.second,
382 "FragmentationAction::performCall() - other atom "
383 +toString(*OtherAtom)+" already present?");
384 }
385 // bonditer follows nicely
386 ASSERT( biniter == position_bins.end(),
387 "FragmentationAction::performCall() - biniter is out of step, it still points at bond "
388 +toString(*biniter)+".");
389 }
390 // and insert
391 globalsaturationpositions.insert(
392 std::make_pair( _atom->getId(),
393 positions_per_neighbor
394 ));
395 }
396 }
397
398 {
399 const enum HydrogenSaturation saturation = params.DoSaturation.get() ? DoSaturate : DontSaturate;
400 const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
401 if (params.types.get().size() != 0) {
402 // store molecule's fragment to file
403 ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
404 exporter.setPrefix(params.prefix.get());
405 exporter.setOutputTypes(params.types.get());
406 exporter();
407 } else {
408 // store molecule's fragment in FragmentJobQueue
409 ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
410 exporter.setLevel(params.level.get());
411 exporter();
412 }
413 }
414
415 // store Adjacency to file
416 {
417 std::string filename = params.prefix.get() + ADJACENCYFILE;
418 std::ofstream AdjacencyFile;
419 AdjacencyFile.open(filename.c_str(), ios::out);
420 AdjacencyList adjacency(atomids);
421 adjacency.StoreToFile(AdjacencyFile);
422 AdjacencyFile.close();
423 }
424
425 World::getInstance().setExitFlag(ExitFlag);
426 end = clock();
427 LOG(0, "STATUS: Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s.");
428
429 return Action::success;
430}
431
432ActionState::ptr FragmentationFragmentationAction::performUndo(ActionState::ptr _state) {
433 return Action::success;
434}
435
436ActionState::ptr FragmentationFragmentationAction::performRedo(ActionState::ptr _state){
437 return Action::success;
438}
439
440bool FragmentationFragmentationAction::canUndo() {
441 return true;
442}
443
444bool FragmentationFragmentationAction::shouldUndo() {
445 return true;
446}
447/** =========== end of function ====================== */
Note: See TracBrowser for help on using the repository browser.