Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/Actions/FragmentationAction/FragmentationAction.cpp

    r5d5550 rbf0ff3  
    4242#include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
    4343#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
    44 #include "Fragmentation/Exporters/SaturatedBond.hpp"
    45 #include "Fragmentation/Exporters/SaturatedFragment.hpp"
    46 #include "Fragmentation/Exporters/SaturationDistanceMaximizer.hpp"
    4744#include "Fragmentation/Fragmentation.hpp"
    4845#include "Fragmentation/Graph.hpp"
     
    9087  // check for selected atoms
    9188  if (world.beginAtomSelection() == world.endAtomSelection()) {
    92     STATUS("There are no atoms selected for fragmentation.");
     89    ELOG(1, "There are no atoms selected for fragmentation.");
    9390    return Action::failure;
    9491  }
     
    264261  }
    265262
    266   // create global saturation positions map
    267   SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
    268   {
    269     // go through each atom
    270     for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
    271         iter != world.endAtomSelection(); ++iter) {
    272       const atom * const _atom = iter->second;
    273 
    274       // skip hydrogens if treated special
    275       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 valence
    282       unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
    283       LOG(3, "DEBUG: There are " << NumberOfPoints
    284           << " places to fill in in total for this atom " << *_atom << ".");
    285 
    286       // check whether there are any bonds with degree larger than 1
    287       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 distances
    297       SaturationDistanceMaximizer::position_bins_t position_bins;
    298       {
    299         // gather all bonds and convert to SaturatedBonds
    300         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 alphas
    310           maximizer();
    311         } else {
    312           // if no higher order bonds, we simply gather the scaled positions
    313         }
    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 map
    319       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               bool
    330               > inserter;
    331           // check whether we treat hydrogen special
    332           if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {
    333             // if hydrogen, forget rescaled position and use original one
    334             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 list
    354           ASSERT (inserter.second,
    355               "FragmentationAction::performCall() - other atom "
    356               +toString(*OtherAtom)+" already present?");
    357         }
    358         // bonditer follows nicely
    359         ASSERT( biniter == position_bins.end(),
    360             "FragmentationAction::performCall() - biniter is out of step, it still points at bond "
    361             +toString(*biniter)+".");
    362       }
    363       // and insert
    364       globalsaturationpositions.insert(
    365           std::make_pair( _atom->getId(),
    366               positions_per_neighbor
    367           ));
    368     }
    369   }
    370 
    371263  {
    372264    const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
     
    374266    if (params.types.get().size() != 0) {
    375267      // store molecule's fragment to file
    376       ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
     268      ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation);
    377269      exporter.setPrefix(params.prefix.get());
    378270      exporter.setOutputTypes(params.types.get());
     
    380272    } else {
    381273      // store molecule's fragment in FragmentJobQueue
    382       ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
     274      ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation);
    383275      exporter.setLevel(params.level.get());
    384276      exporter();
Note: See TracChangeset for help on using the changeset viewer.