- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Actions/FragmentationAction/FragmentationAction.cpp
r5d5550 rbf0ff3 42 42 #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp" 43 43 #include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp" 44 #include "Fragmentation/Exporters/SaturatedBond.hpp"45 #include "Fragmentation/Exporters/SaturatedFragment.hpp"46 #include "Fragmentation/Exporters/SaturationDistanceMaximizer.hpp"47 44 #include "Fragmentation/Fragmentation.hpp" 48 45 #include "Fragmentation/Graph.hpp" … … 90 87 // check for selected atoms 91 88 if (world.beginAtomSelection() == world.endAtomSelection()) { 92 STATUS("There are no atoms selected for fragmentation.");89 ELOG(1, "There are no atoms selected for fragmentation."); 93 90 return Action::failure; 94 91 } … … 264 261 } 265 262 266 // create global saturation positions map267 SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;268 {269 // go through each atom270 for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();271 iter != world.endAtomSelection(); ++iter) {272 const atom * const _atom = iter->second;273 274 // skip hydrogens if treated special275 const enum HydrogenTreatment treatment = params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;276 if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {277 LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);278 continue;279 }280 281 // get the valence282 unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();283 LOG(3, "DEBUG: There are " << NumberOfPoints284 << " places to fill in in total for this atom " << *_atom << ".");285 286 // check whether there are any bonds with degree larger than 1287 unsigned int SumOfDegrees = 0;288 bool PresentHigherBonds = false;289 const BondList &bondlist = _atom->getListOfBonds();290 for (BondList::const_iterator bonditer = bondlist.begin();291 bonditer != bondlist.end(); ++bonditer) {292 SumOfDegrees += (*bonditer)->getDegree();293 PresentHigherBonds |= (*bonditer)->getDegree() > 1;294 }295 296 // check whether there are alphas to maximize the hydrogens distances297 SaturationDistanceMaximizer::position_bins_t position_bins;298 {299 // gather all bonds and convert to SaturatedBonds300 SaturationDistanceMaximizer::PositionContainers_t CutBonds;301 for (BondList::const_iterator bonditer = bondlist.begin();302 bonditer != bondlist.end(); ++bonditer) {303 CutBonds.push_back(304 SaturatedBond::ptr(new SaturatedBond(*(bonditer->get()), *_atom) )305 );306 }307 SaturationDistanceMaximizer maximizer(CutBonds);308 if (PresentHigherBonds) {309 // then find best alphas310 maximizer();311 } else {312 // if no higher order bonds, we simply gather the scaled positions313 }314 position_bins = maximizer.getAllPositionBins();315 LOG(4, "DEBUG: Positions for atom " << *_atom << " are " << position_bins);316 }317 318 // convert into the desired entry in the map319 SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;320 {321 BondList::const_iterator bonditer = bondlist.begin();322 SaturationDistanceMaximizer::position_bins_t::const_iterator biniter =323 position_bins.begin();324 325 for (;bonditer != bondlist.end(); ++bonditer, ++biniter) {326 const atom * const OtherAtom = (*bonditer)->GetOtherAtom(_atom);327 std::pair<328 SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,329 bool330 > inserter;331 // check whether we treat hydrogen special332 if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {333 // if hydrogen, forget rescaled position and use original one334 inserter =335 positions_per_neighbor.insert(336 std::make_pair(337 OtherAtom->getId(),338 SaturatedFragment::SaturationsPositions_t(339 1, OtherAtom->getPosition() - _atom->getPosition())340 )341 );342 } else {343 inserter =344 positions_per_neighbor.insert(345 std::make_pair(346 OtherAtom->getId(),347 SaturatedFragment::SaturationsPositions_t(348 biniter->begin(),349 biniter->end())350 )351 );352 }353 // if already pressent, add to this present list354 ASSERT (inserter.second,355 "FragmentationAction::performCall() - other atom "356 +toString(*OtherAtom)+" already present?");357 }358 // bonditer follows nicely359 ASSERT( biniter == position_bins.end(),360 "FragmentationAction::performCall() - biniter is out of step, it still points at bond "361 +toString(*biniter)+".");362 }363 // and insert364 globalsaturationpositions.insert(365 std::make_pair( _atom->getId(),366 positions_per_neighbor367 ));368 }369 }370 371 263 { 372 264 const enum HydrogenSaturation saturation = params.DoSaturation.get() ? DoSaturate : DontSaturate; … … 374 266 if (params.types.get().size() != 0) { 375 267 // store molecule's fragment to file 376 ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation , globalsaturationpositions);268 ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation); 377 269 exporter.setPrefix(params.prefix.get()); 378 270 exporter.setOutputTypes(params.types.get()); … … 380 272 } else { 381 273 // store molecule's fragment in FragmentJobQueue 382 ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation , globalsaturationpositions);274 ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation); 383 275 exporter.setLevel(params.level.get()); 384 276 exporter();
Note:
See TracChangeset
for help on using the changeset viewer.