source: src/Fragmentation/Fragmentation.cpp@ 9511c7

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * Fragmentation.cpp
 *
 *  Created on: Oct 18, 2011
 *      Author: heber
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "Fragmentation.hpp"

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Descriptors/MoleculeDescriptor.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/AdaptivityMap.hpp"
#include "Fragmentation/AtomMask.hpp"
#include "Fragmentation/fragmentation_helpers.hpp"
#include "Fragmentation/Graph.hpp"
#include "Fragmentation/helpers.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/PowerSetGenerator.hpp"
#include "Fragmentation/UniqueFragments.hpp"
#include "Graph/BondGraph.hpp"
#include "Graph/CheckAgainstAdjacencyFile.hpp"
#include "Graph/ListOfLocalAtoms.hpp"
#include "molecule.hpp"
#include "MoleculeLeafClass.hpp"
#include "MoleculeListClass.hpp"
#include "Parser/FormatParserStorage.hpp"
#include "World.hpp"


/** Constructor of class Fragmentation.
 *
 * \param _mol molecule for internal use (looking up atoms)
 * \param _FileChecker instance contains adjacency parsed from elsewhere
 * \param _saturation whether to treat hydrogen special and saturate dangling bonds or not
 */
Fragmentation::Fragmentation(molecule *_mol, CheckAgainstAdjacencyFile &_FileChecker, const enum HydrogenSaturation _saturation) :
  mol(_mol),
  saturation(_saturation),
  FileChecker(_FileChecker)
{}

/** Destructor of class Fragmentation.
 *
 */
Fragmentation::~Fragmentation()
{}


/** Performs a many-body bond order analysis for a given bond order.
 * -# parses adjacency, keysets and orderatsite files
 * -# RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energ
y contribution", and that's why this consciously not done in the following loop)
 * -# in a loop over all subgraphs
 *  -# calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
 *  -# creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)
 * -# combines the generated molecule lists from all subgraphs
 * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
 * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
 * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
 * subgraph in the MoleculeListClass.
 * \param atomids atoms to fragment, used in AtomMask
 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
 * \param prefix prefix string for every fragment file name (may include path)
 * \param &DFS contains bond structure analysis data
 * \return 1 - continue, 2 - stop (no fragmentation occured)
 */
int Fragmentation::FragmentMolecule(const std::vector<atomId_t> &atomids, int Order, std::string prefix, DepthFirstSearchAnalysis &DFS)
{
  MoleculeListClass *BondFragments = NULL;
  int FragmentCounter;
  MoleculeLeafClass *MolecularWalker = NULL;
  MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
  fstream File;
  bool FragmentationToDo = true;
  bool CheckOrder = false;
  Graph **FragmentList = NULL;
  Graph TotalGraph;     // graph with all keysets however local numbers
  int TotalNumberOfKeySets = 0;
  AtomMask_t AtomMask(atomids);

  LOG(0, endl);
  switch (saturation) {
    case DoSaturate:
      LOG(0, "I will treat hydrogen special and saturate dangling bonds with it.");
      break;
    case DontSaturate:
      LOG(0, "Hydrogen is treated just like the rest of the lot.");
      break;
    default:
      ASSERT(0, "Fragmentation::FragmentMolecule() - there is a saturation setting which I have no idea about.");
      break;
  }

  // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++

  // ===== 1. Check whether bond structure is same as stored in files ====

  // === compare it with adjacency file ===
  {
    const std::vector<atomId_t> globalids = getGlobalIdsFromLocalIds(*mol, atomids);
    FragmentationToDo = FragmentationToDo && FileChecker(globalids);
  }

  // === reset bond degree and perform CorrectBondDegree ===
  for(World::MoleculeIterator iter = World::getInstance().getMoleculeIter();
      iter != World::getInstance().moleculeEnd();
      ++iter) {
    // correct bond degree
    World::AtomComposite Set = (*iter)->getAtomSet();
    World::getInstance().getBondGraph()->CorrectBondDegree(Set);
  }

  // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
  // NOTE: We assume here that DFS has been performed and molecule structure is current.
  Subgraphs = DFS.getMoleculeStructure();

  // ===== 3. if structure still valid, parse key set file and others =====
  Graph ParsedFragmentList;
  FragmentationToDo = FragmentationToDo && ParsedFragmentList.ParseKeySetFile(prefix);

  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(atomids, prefix);

  // =================================== Begin of FRAGMENTATION ===============================
  // ===== 6a. assign each keyset to its respective subgraph =====
  const int MolCount = World::getInstance().numMolecules();
  FragmentCounter = 0;
  {
    ListOfLocalAtoms_t *ListOfLocalAtoms = new ListOfLocalAtoms_t[MolCount];
    Subgraphs->next->AssignKeySetsToFragment(mol, &ParsedFragmentList, ListOfLocalAtoms, FragmentList, FragmentCounter, true);
    delete[] ListOfLocalAtoms;
  }

  // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
  bool LoopDoneAlready = false;
  while ((CheckOrder = CheckOrderAtSite(AtomMask, &ParsedFragmentList, Order, prefix, LoopDoneAlready))) {
    FragmentationToDo = FragmentationToDo || CheckOrder;
    LoopDoneAlready = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0), saturation);

    // ===== 7. fill the bond fragment list =====
    FragmentCounter = 0;
    MolecularWalker = Subgraphs;
    while (MolecularWalker->next != NULL) {
      MolecularWalker = MolecularWalker->next;
      LOG(1, "Fragmenting subgraph " << MolecularWalker << ".");
      if (MolecularWalker->Leaf->hasBondStructure()) {
        // call BOSSANOVA method
        LOG(0, endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= ");
        FragmentBOSSANOVA(MolecularWalker->Leaf, FragmentList[FragmentCounter], RootStack[FragmentCounter]);
      } else {
        ELOG(1, "Subgraph " << MolecularWalker << " has no atoms!");
      }
      FragmentCounter++;  // next fragment list
    }
  }
  LOG(2, "CheckOrder is " << CheckOrder << ".");
  delete[](RootStack);

  // ==================================== End of FRAGMENTATION ============================================

  // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
  Subgraphs->next->TranslateIndicesToGlobalIDs(FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);

  // free subgraph memory again
  FragmentCounter = 0;
  while (Subgraphs != NULL) {
    // remove entry in fragment list
    // remove subgraph fragment
    MolecularWalker = Subgraphs->next;
    Subgraphs->Leaf = NULL;
    delete(Subgraphs);
    Subgraphs = MolecularWalker;
  }
  // free fragment list
  for (int i=0; i< FragmentCounter; ++i )
    delete(FragmentList[i]);
  delete[](FragmentList);

  LOG(0, FragmentCounter << " subgraph fragments have been removed.");

  // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
  //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
  // allocate memory for the pointer array and transmorph graphs into full molecular fragments
  BondFragments = new MoleculeListClass(World::getPointer());
  int k=0;
  for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
    KeySet test = (*runner).first;
    LOG(0, "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << ".");
    BondFragments->insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
    k++;
  }
  LOG(0, k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets.");

  // ===== 9. Save fragments' configuration and keyset files et al to disk ===
  if (BondFragments->ListOfMolecules.size() != 0) {
    // create the SortIndex from BFS labels to order in the config file
    std::map<atomId_t, int> SortIndex;
    CreateMappingLabelsToConfigSequence(SortIndex);

    LOG(1, "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs");
    bool write_status = true;
    for (std::vector<std::string>::const_iterator iter = typelist.begin();
        iter != typelist.end();
        ++iter) {
      LOG(2, "INFO: Writing bond fragments for type " << (*iter) << ".");
      write_status = write_status
      && BondFragments->OutputConfigForListOfFragments(
          prefix,
          FormatParserStorage::getInstance().getTypeFromName(*iter));
    }
    if (write_status)
      LOG(1, "All configs written.");
    else
      LOG(1, "Some config writing failed.");

    // store force index reference file
    BondFragments->StoreForcesFile(prefix, SortIndex);

    // store keysets file
    TotalGraph.StoreKeySetFile(prefix);

    {
      // store Adjacency file
      std::string filename = prefix + ADJACENCYFILE;
      mol->StoreAdjacencyToFile(filename);
    }

    // store Hydrogen saturation correction file
    BondFragments->AddHydrogenCorrection(prefix);

    // store adaptive orders into file
    StoreOrderAtSiteFile(prefix);

    // restore orbital and Stop values
    //CalculateOrbitals(*configuration);
  } else {
    LOG(1, "FragmentList is zero on return, splitting failed.");
  }
  // remove all create molecules again from the World including their atoms
  for (MoleculeList::iterator iter = BondFragments->ListOfMolecules.begin();
      !BondFragments->ListOfMolecules.empty();
      iter = BondFragments->ListOfMolecules.begin()) {
    // remove copied atoms and molecule again
    molecule *mol = *iter;
    mol->removeAtomsinMolecule();
    World::getInstance().destroyMolecule(mol);
    BondFragments->ListOfMolecules.erase(iter);
  }
  delete(BondFragments);
  LOG(0, "End of bond fragmentation.");

  return ((int)(!FragmentationToDo)+1);    // 1 - continue, 2 - stop (no fragmentation occured)
};


/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
 * -# constructs a complete keyset of the molecule
 * -# In a loop over all possible roots from the given rootstack
 *  -# increases order of root site
 *  -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
 *  -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
as the restricted one and each site in the set as the root)
 *  -# these are merged into a fragment list of keysets
 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
 * Important only is that we create all fragments, it is not important if we create them more than once
 * as these copies are filtered out via use of the hash table (KeySet).
 * \param *out output stream for debugging
 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
 * \return pointer to Graph list
 */
void Fragmentation::FragmentBOSSANOVA(molecule *mol, Graph *&FragmentList, KeyStack &RootStack)
{
  Graph ***FragmentLowerOrdersList = NULL;
  int NumLevels = 0;
  int NumMolecules = 0;
  int TotalNumMolecules = 0;
  int *NumMoleculesOfOrder = NULL;
  int Order = 0;
  int UpgradeCount = RootStack.size();
  KeyStack FragmentRootStack;
  int RootKeyNr = 0;
  int RootNr = 0;
  UniqueFragments FragmentSearch;

  LOG(0, "Begin of FragmentBOSSANOVA.");

  // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
  // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
  NumMoleculesOfOrder = new int[UpgradeCount];
  FragmentLowerOrdersList = new Graph**[UpgradeCount];

  for(int i=0;i<UpgradeCount;i++) {
    NumMoleculesOfOrder[i] = 0;
    FragmentLowerOrdersList[i] = NULL;
  }

  FragmentSearch.Init(mol->FindAtom(RootKeyNr));

  // Construct the complete KeySet which we need for topmost level only (but for all Roots)
  KeySet CompleteMolecule;
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    CompleteMolecule.insert((*iter)->GetTrueFather()->getNr());
  }

  // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
  // each has to be split up. E.g. for the second level we have one from 5th, one from 4th, two from 3th (which in turn is one from 5th, one from 4th),
  // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
  // with the order along the cells as this: 5433222211111111 for BondOrder 5 needing 16=pow(2,5-1) cells (only we use bit-shifting which is faster)
  RootNr = 0;   // counts through the roots in RootStack
  while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
    RootKeyNr = RootStack.front();
    RootStack.pop_front();
    atom *Walker = mol->FindAtom(RootKeyNr);
    // check cyclic lengths
    //if ((MinimumRingSize[Walker->GetTrueFather()->getNr()] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->getNr()])) {
    //  LOG(0, "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed.");
    //} else
    {
      // increase adaptive order by one
      Walker->GetTrueFather()->AdaptiveOrder++;
      Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;

      // initialise Order-dependent entries of UniqueFragments structure
      class PowerSetGenerator PSG(&FragmentSearch, Walker->AdaptiveOrder);

      // allocate memory for all lower level orders in this 1D-array of ptrs
      NumLevels = 1 << (Order-1); // (int)pow(2,Order);
      FragmentLowerOrdersList[RootNr] = new Graph*[NumLevels];
      for (int i=0;i<NumLevels;i++)
        FragmentLowerOrdersList[RootNr][i] = NULL;

      // create top order where nothing is reduced
      LOG(0, "==============================================================================================================");
      LOG(0, "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining."); // , NumLevels is " << NumLevels << "

      // Create list of Graphs of current Bond Order (i.e. F_{ij})
      FragmentLowerOrdersList[RootNr][0] = new Graph;
      FragmentSearch.PrepareForPowersetGeneration(1., FragmentLowerOrdersList[RootNr][0], Walker);
      NumMoleculesOfOrder[RootNr] = PSG(CompleteMolecule, saturation);

      // output resulting number
      LOG(1, "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << ".");
      if (NumMoleculesOfOrder[RootNr] != 0) {
        NumMolecules = 0;
      } else {
        Walker->GetTrueFather()->MaxOrder = true;
      }
      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
      //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
      TotalNumMolecules += NumMoleculesOfOrder[RootNr];
//      LOG(1, "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << ".");
      RootStack.push_back(RootKeyNr); // put back on stack
      RootNr++;
    }
  }
  LOG(0, "==============================================================================================================");
  LOG(1, "Total number of resulting molecules is: " << TotalNumMolecules << ".");
  LOG(0, "==============================================================================================================");

  // cleanup FragmentSearch structure
  FragmentSearch.Cleanup();

  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
  // 5433222211111111
  // 43221111
  // 3211
  // 21
  // 1

  // Subsequently, we combine all into a single list (FragmentList)
  CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, mol);
  FreeAllOrdersList(FragmentLowerOrdersList, RootStack, mol);
  delete[](NumMoleculesOfOrder);

  LOG(0, "End of FragmentBOSSANOVA.");
};

/** Stores a fragment from \a KeySet into \a molecule.
 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
 * molecule and adds missing hydrogen where bonds were cut.
 * \param *out output stream for debugging messages
 * \param &Leaflet pointer to KeySet structure
 * \param IsAngstroem whether we have Ansgtroem or bohrradius
 * \return pointer to constructed molecule
 */
molecule * Fragmentation::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
{
  Info info(__func__);
  ListOfLocalAtoms_t SonList;
  molecule *Leaf = World::getInstance().createMolecule();

  StoreFragmentFromKeySet_Init(mol, Leaf, Leaflet, SonList);
  // create the bonds between all: Make it an induced subgraph and add hydrogen
//  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
  CreateInducedSubgraphOfFragment(mol, Leaf, SonList, IsAngstroem);

  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
  return Leaf;
};


/** Estimates by educated guessing (using upper limit) the expected number of fragments.
 * The upper limit is
 * \f[
 *  n = N \cdot C^k
 * \f]
 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
 * \param *out output stream for debugging
 * \param order bond order k
 * \return number n of fragments
 */
int Fragmentation::GuesstimateFragmentCount(int order)
{
  size_t c = 0;
  int FragmentCount;
  // get maximum bond degree
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    const BondList& ListOfBonds = (*iter)->getListOfBonds();
    c = (ListOfBonds.size() > c) ? ListOfBonds.size() : c;
  }
  FragmentCount = (saturation == DoSaturate ? mol->getNoNonHydrogen() : mol->getAtomCount()) *(1 << (c*order));
  LOG(1, "Upper limit for this subgraph is " << FragmentCount << " for "
      << mol->getNoNonHydrogen() << " non-H atoms with maximum bond degree of " << c << ".");
  return FragmentCount;
};


/** Checks whether the OrderAtSite is still below \a Order at some site.
 * \param AtomMask defines true/false per global Atom::Nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
 * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
 * \param LoopDoneAlready indicate whether we have done a fragmentation loop already
 * \return true - needs further fragmentation, false - does not need fragmentation
 */
bool Fragmentation::CheckOrderAtSite(AtomMask_t &AtomMask, Graph *GlobalKeySetList, int Order, std::string path, bool LoopDoneAlready)
{
  bool status = false;

  // initialize mask list
  AtomMask.clear();

  if (Order < 0) { // adaptive increase of BondOrder per site
    if (LoopDoneAlready)  // break after one step
      return false;

    // transmorph graph keyset list into indexed KeySetList
    if (GlobalKeySetList == NULL) {
      ELOG(1, "Given global key set list (graph) is NULL!");
      return false;
    }
    AdaptivityMap * IndexKeySetList = GlobalKeySetList->GraphToAdaptivityMap();

    // parse the EnergyPerFragment file
    IndexKeySetList->ScanAdaptiveFileIntoMap(path); // (Root No., (Value, Order)) !
    // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
    IndexKeySetList->ReMapAdaptiveCriteriaListToValue(mol);

    // pick the ones still below threshold and mark as to be adaptively updated
    if (IndexKeySetList->IsAdaptiveCriteriaListEmpty()) {
      ELOG(2, "Unable to parse file, incrementing all.");
      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
        if ((saturation == DontSaturate) || ((*iter)->getType()->getAtomicNumber() != 1)) // skip hydrogen
        {
          AtomMask.setTrue((*iter)->getNr());  // include all (non-hydrogen) atoms
          status = true;
        }
      }
    } else {
      IndexKeySetList->MarkUpdateCandidates(AtomMask, Order, mol);
    }

    delete[](IndexKeySetList);
  } else { // global increase of Bond Order
    for(molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
      if ((saturation == DontSaturate) || ((*iter)->getType()->getAtomicNumber() != 1)) // skip hydrogen
      {
        AtomMask.setTrue((*iter)->getNr());  // include all (non-hydrogen) atoms
        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->getNr()]))
          status = true;
      }
    }
    if ((!Order) && (!LoopDoneAlready))  // single stepping, just check
      status = true;

    if (!status) {
      if (Order == 0)
        LOG(1, "Single stepping done.");
      else
        LOG(1, "Order at every site is already equal or above desired order " << Order << ".");
    }
  }

  PrintAtomMask(AtomMask, mol->getAtomCount()); // for debugging

  return status;
};

/** Stores pairs (Atom::Nr, Atom::AdaptiveOrder) into file.
 * Atoms not present in the file get "-1".
 * \param &path path to file ORDERATSITEFILE
 * \return true - file writable, false - not writable
 */
bool Fragmentation::StoreOrderAtSiteFile(std::string &path)
{
  string line;
  ofstream file;

  line = path + ORDERATSITEFILE;
  file.open(line.c_str());
  LOG(1, "Writing OrderAtSite " << ORDERATSITEFILE << " ... ");
  if (file.good()) {
    for_each(mol->begin(),mol->end(),bind2nd(mem_fun(&atom::OutputOrder), &file));
    file.close();
    LOG(1, "done.");
    return true;
  } else {
    LOG(1, "failed to open file " << line << ".");
    return false;
  }
};


/** Parses pairs(Atom::Nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
 * Atoms not present in the file get "0".
 * \param atomids atoms to fragment, used in AtomMask
 * \param &path path to file ORDERATSITEFILEe
 * \return true - file found and scanned, false - file not found
 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
 */
bool Fragmentation::ParseOrderAtSiteFromFile(const std::vector<atomId_t> &atomids, std::string &path)
{
  typedef unsigned char order_t;
  typedef std::map<atomId_t, order_t> OrderArray_t;
  OrderArray_t OrderArray;
  AtomMask_t MaxArray(atomids);
  bool status;
  int AtomNr, value;
  string line;
  ifstream file;

  LOG(1, "Begin of ParseOrderAtSiteFromFile");
  line = path + ORDERATSITEFILE;
  file.open(line.c_str());
  if (file.good()) {
    while (!file.eof()) { // parse from file
      AtomNr = -1;
      file >> AtomNr;
      if (AtomNr != -1) {   // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
        file >> value;
        OrderArray[AtomNr] = value;
        file >> value;
        MaxArray.setValue(AtomNr, (bool)value);
        //LOG(2, "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << ".");
      }
    }
    file.close();

    // set atom values
    for(molecule::iterator iter=mol->begin();iter!=mol->end();++iter){
      (*iter)->AdaptiveOrder = OrderArray[(*iter)->getNr()];
      (*iter)->MaxOrder = MaxArray.isTrue((*iter)->getNr());
    }
    //SetAtomValueToIndexedArray( OrderArray, &atom::getNr(), &atom::AdaptiveOrder );
    //SetAtomValueToIndexedArray( MaxArray, &atom::getNr(), &atom::MaxOrder );

    LOG(1, "\t ... done.");
    status = true;
  } else {
    LOG(1, "\t ... failed to open file " << line << ".");
    status = false;
  }

  LOG(1, "End of ParseOrderAtSiteFromFile");
  return status;
};

/** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
 * \param &SortIndex Mapping array of size molecule::AtomCount
 * \return true - success, false - failure of SortIndex alloc
 */
bool Fragmentation::CreateMappingLabelsToConfigSequence(std::map<atomId_t, int> &SortIndex)
{
  if (!SortIndex.empty()) {
    LOG(1, "SortIndex has " << SortIndex.size() << " entries and is not empty as expected.");
    return false;
  }

  int AtomNo = 0;
  for(molecule::const_iterator iter=mol->begin();iter!=mol->end();++iter){
    const int id = (*iter)->getNr();
#ifndef NDEBUG
    std::pair<std::map<atomId_t, int>::const_iterator, bool> inserter =
#endif
    SortIndex.insert( std::make_pair(id, AtomNo++) );
    ASSERT( inserter.second ,
        "Fragmentation::CreateMappingLabelsToConfigSequence() - same SortIndex set twice.");
  }

  return true;
};


/** Initializes some value for putting fragment of \a *mol into \a *Leaf.
 * \param *mol total molecule
 * \param *Leaf fragment molecule
 * \param &Leaflet pointer to KeySet structure
 * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
 * \return number of atoms in fragment
 */
int Fragmentation::StoreFragmentFromKeySet_Init(molecule *mol, molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
{
  atom *FatherOfRunner = NULL;

  // first create the minimal set of atoms from the KeySet
  int size = 0;
  for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
    FatherOfRunner = mol->FindAtom((*runner));  // find the id
    SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
    size++;
  }
  return size;
};


/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
 * \param *out output stream for debugging messages
 * \param *mol total molecule
 * \param *Leaf fragment molecule
 * \param IsAngstroem whether we have Ansgtroem or bohrradius
 * \param SonList list which atom of \a *Leaf is a son of which atom in \a *mol
 */
void Fragmentation::CreateInducedSubgraphOfFragment(molecule *mol, molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
{
  bool LonelyFlag = false;
  atom *OtherFather = NULL;
  atom *FatherOfRunner = NULL;

  // we increment the iter just before skipping the hydrogen
  // as we use AddBond, we cannot have a const_iterator here
  for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
    LonelyFlag = true;
    FatherOfRunner = (*iter)->father;
    ASSERT(FatherOfRunner,"Atom without father found");
    if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
      // create all bonds
      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
          BondRunner != ListOfBonds.end();
          ++BondRunner) {
        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
        if (SonList.find(OtherFather->getNr()) != SonList.end()) {
//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
//              << " is bound to " << *OtherFather << ", whose son is "
//              << *SonList[OtherFather->getNr()] << ".");
          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
            std::stringstream output;
//            output << "ACCEPT: Adding Bond: "
            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->BondDegree);
//            LOG(3, output.str());
            //NumBonds[(*iter)->getNr()]++;
          } else {
//            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
          }
          LonelyFlag = false;
        } else {
//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
//              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
          if (saturation == DoSaturate) {
//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
              exit(1);
          }
          //NumBonds[(*iter)->getNr()] += Binder->BondDegree;
        }
      }
    } else {
    ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
    }
    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
      LOG(0, **iter << "has got bonds only to hydrogens!");
    }
    ++iter;
    if (saturation == DoSaturate) {
      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
        iter++;
      }
    }
  }
};

/** Sets the desired output types of the fragment configurations.
 *
 * @param types vector of desired types.
 */
void Fragmentation::setOutputTypes(const std::vector<std::string> &types)
{
  typelist = types;
}