Changeset 246e13 for src/Fragmentation


Ignore:
Timestamp:
Oct 25, 2011, 12:08:03 PM (13 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, 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_BoundInBox_CenterInBox_MoleculeActions, 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, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, 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
Children:
f0674a
Parents:
212c179
git-author:
Frederik Heber <heber@…> (10/19/11 17:07:28)
git-committer:
Frederik Heber <heber@…> (10/25/11 12:08:03)
Message:

Placed FragmentMolecule, FragmentBOSSANOVA and subfunctions into own class Fragmentation.

Location:
src/Fragmentation
Files:
1 added
1 edited
1 moved

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Fragmentation.cpp

    r212c179 r246e13  
    77
    88/*
    9  * molecule_fragmentation.cpp
     9 * Fragmentation.cpp
    1010 *
    11  *  Created on: Oct 5, 2009
     11 *  Created on: Oct 18, 2011
    1212 *      Author: heber
    1313 */
    1414
    15 // include config.h
    1615#ifdef HAVE_CONFIG_H
    1716#include <config.h>
    … …  
    2019#include "CodePatterns/MemDebug.hpp"
    2120
    22 #include <cstring>
     21#include "Fragmentation.hpp"
     22
     23#include "CodePatterns/Assert.hpp"
     24#include "CodePatterns/Log.hpp"
    2325
    2426#include "atom.hpp"
    2527#include "Bond/bond.hpp"
    26 #include "Box.hpp"
    27 #include "CodePatterns/Verbose.hpp"
    28 #include "CodePatterns/Log.hpp"
    29 #include "config.hpp"
    3028#include "Element/element.hpp"
    3129#include "Element/periodentafel.hpp"
    … …  
    3331#include "Fragmentation/PowerSetGenerator.hpp"
    3432#include "Fragmentation/UniqueFragments.hpp"
     33#include "graph.hpp"
    3534#include "Graph/BondGraph.hpp"
    3635#include "Graph/CheckAgainstAdjacencyFile.hpp"
    37 #include "Graph/CyclicStructureAnalysis.hpp"
    38 #include "Graph/DepthFirstSearchAnalysis.hpp"
    39 #include "Helpers/helpers.hpp"
    40 #include "LinearAlgebra/RealSpaceMatrix.hpp"
    4136#include "molecule.hpp"
    4237#include "World.hpp"
    4338
    44 /************************************* Functions for class molecule *********************************/
    45 
    46 
    47 /** Estimates by educated guessing (using upper limit) the expected number of fragments.
    48  * The upper limit is
    49  * \f[
    50  *  n = N \cdot C^k
    51  * \f]
    52  * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
    53  * \param *out output stream for debugging
    54  * \param order bond order k
    55  * \return number n of fragments
    56  */
    57 int molecule::GuesstimateFragmentCount(int order)
    58 {
    59   size_t c = 0;
    60   int FragmentCount;
    61   // get maximum bond degree
    62   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    63     const BondList& ListOfBonds = (*iter)->getListOfBonds();
    64     c = (ListOfBonds.size() > c) ? ListOfBonds.size() : c;
    65   }
    66   FragmentCount = NoNonHydrogen*(1 << (c*order));
    67   DoLog(1) && (Log() << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl);
    68   return FragmentCount;
    69 };
    70 
    71 /** Checks whether the OrderAtSite is still below \a Order at some site.
    72  * \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
    73  * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
    74  * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
    75  * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
    76  * \return true - needs further fragmentation, false - does not need fragmentation
    77  */
    78 bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, std::string path)
    79 {
    80   bool status = false;
    81 
    82   // initialize mask list
    83   for(int i=getAtomCount();i--;)
    84     AtomMask[i] = false;
    85 
    86   if (Order < 0) { // adaptive increase of BondOrder per site
    87     if (AtomMask[getAtomCount()] == true)  // break after one step
    88       return false;
    89 
    90     // transmorph graph keyset list into indexed KeySetList
    91     if (GlobalKeySetList == NULL) {
    92       DoeLog(1) && (eLog()<< Verbose(1) << "Given global key set list (graph) is NULL!" << endl);
    93       return false;
    94     }
    95     map<int,KeySet> *IndexKeySetList = GraphToIndexedKeySet(GlobalKeySetList);
    96 
    97     // parse the EnergyPerFragment file
    98     map<int, pair<double,int> > *AdaptiveCriteriaList = ScanAdaptiveFileIntoMap(path, *IndexKeySetList); // (Root No., (Value, Order)) !
    99     if (AdaptiveCriteriaList->empty()) {
    100       DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
    101       for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    102     #ifdef ADDHYDROGEN
    103         if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
    104     #endif
    105         {
    106           AtomMask[(*iter)->getNr()] = true;  // include all (non-hydrogen) atoms
    107           status = true;
    108         }
    109       }
    110     }
    111     // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
    112     map<double, pair<int,int> > *FinalRootCandidates = ReMapAdaptiveCriteriaListToValue(AdaptiveCriteriaList, this);
    113 
    114     // pick the ones still below threshold and mark as to be adaptively updated
    115     MarkUpdateCandidates(AtomMask, *FinalRootCandidates, Order, this);
    116 
    117     delete[](IndexKeySetList);
    118     delete[](AdaptiveCriteriaList);
    119     delete[](FinalRootCandidates);
    120   } else { // global increase of Bond Order
    121     for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    122   #ifdef ADDHYDROGEN
    123       if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
    124   #endif
    125       {
    126         AtomMask[(*iter)->getNr()] = true;  // include all (non-hydrogen) atoms
    127         if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->getNr()]))
    128           status = true;
    129       }
    130     }
    131     if ((!Order) && (!AtomMask[getAtomCount()]))  // single stepping, just check
    132       status = true;
    133 
    134     if (!status) {
    135       if (Order == 0)
    136         DoLog(1) && (Log() << Verbose(1) << "Single stepping done." << endl);
    137       else
    138         DoLog(1) && (Log() << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl);
    139     }
    140   }
    141 
    142   PrintAtomMask(AtomMask, getAtomCount()); // for debugging
    143 
    144   return status;
    145 };
    146 
    147 /** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
    148  * \param *out output stream for debugging
    149  * \param *&SortIndex Mapping array of size molecule::AtomCount
    150  * \return true - success, false - failure of SortIndex alloc
    151  */
    152 bool molecule::CreateMappingLabelsToConfigSequence(int *&SortIndex)
    153 {
    154   if (SortIndex != NULL) {
    155     DoLog(1) && (Log() << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl);
    156     return false;
    157   }
    158   SortIndex = new int[getAtomCount()];
    159   for(int i=getAtomCount();i--;)
    160     SortIndex[i] = -1;
    161 
    162   int AtomNo = 0;
    163   for(internal_iterator iter=atoms.begin();iter!=atoms.end();++iter){
    164     ASSERT(SortIndex[(*iter)->getNr()]==-1,"Same SortIndex set twice");
    165     SortIndex[(*iter)->getNr()] = AtomNo++;
    166   }
    167 
    168   return true;
    169 };
    170 
    171 
    172 
    173 /** Creates a lookup table for true father's Atom::Nr -> atom ptr.
    174  * \param *start begin of list (STL iterator, i.e. first item)
    175  * \paran *end end of list (STL iterator, i.e. one past last item)
    176  * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
    177  * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
    178  * \return true - success, false - failure
    179  */
    180 bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
    181 {
    182   bool status = true;
    183   int AtomNo;
    184 
    185   if (LookupTable != NULL) {
    186     Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
    187     return false;
    188   }
    189 
    190   // count them
    191   if (count == 0) {
    192     for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::Nr -> atom) used as a marker table lateron
    193       count = (count < (*iter)->GetTrueFather()->getNr()) ? (*iter)->GetTrueFather()->getNr() : count;
    194     }
    195   }
    196   if (count <= 0) {
    197     Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
    198     return false;
    199   }
    200 
    201   // allocate and fill
    202   LookupTable = new atom *[count];
    203   if (LookupTable == NULL) {
    204     eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
    205     performCriticalExit();
    206     status = false;
    207   } else {
    208     for (int i=0;i<count;i++)
    209       LookupTable[i] = NULL;
    210     for (molecule::iterator iter = begin(); iter != end(); ++iter) {
    211       AtomNo = (*iter)->GetTrueFather()->getNr();
    212       if ((AtomNo >= 0) && (AtomNo < count)) {
    213         //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl;
    214         LookupTable[AtomNo] = (*iter);
    215       } else {
    216         Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl;
    217         status = false;
    218         break;
    219       }
    220     }
    221   }
    222 
    223   return status;
    224 };
     39
     40/** Constructor of class Fragmentation.
     41 *
     42 */
     43Fragmentation::Fragmentation(molecule *_mol) :
     44  mol(_mol)
     45{}
     46
     47/** Destructor of class Fragmentation.
     48 *
     49 */
     50Fragmentation::~Fragmentation()
     51{}
    22552
    22653
    … …  
    24370 * \return 1 - continue, 2 - stop (no fragmentation occured)
    24471 */
    245 int molecule::FragmentMolecule(int Order, std::string &prefix, DepthFirstSearchAnalysis &DFS)
     72int Fragmentation::FragmentMolecule(int Order, std::string &prefix, DepthFirstSearchAnalysis &DFS)
    24673{
    24774  MoleculeListClass *BondFragments = NULL;
    … …  
    308135    ListOfLocalAtoms[i] = NULL;
    309136  FragmentCounter = 0;
    310   Subgraphs->next->AssignKeySetsToFragment(this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, FragmentCounter, true);
     137  Subgraphs->next->AssignKeySetsToFragment(mol, ParsedFragmentList, ListOfLocalAtoms, FragmentList, FragmentCounter, true);
    311138  delete[](ListOfLocalAtoms);
    312139
    313140  // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
    314141  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
    315   AtomMask = new bool[getAtomCount()+1];
    316   AtomMask[getAtomCount()] = false;
     142  AtomMask = new bool[mol->getAtomCount()+1];
     143  AtomMask[mol->getAtomCount()] = false;
    317144  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
    318145  while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, prefix))) {
    319146    FragmentationToDo = FragmentationToDo || CheckOrder;
    320     AtomMask[getAtomCount()] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
     147    AtomMask[mol->getAtomCount()] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    321148    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    322149    Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0));
    … …  
    331158        // call BOSSANOVA method
    332159        DoLog(0) && (Log() << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl);
    333         MolecularWalker->Leaf->FragmentBOSSANOVA(FragmentList[FragmentCounter], RootStack[FragmentCounter]);
     160        FragmentBOSSANOVA(MolecularWalker->Leaf, FragmentList[FragmentCounter], RootStack[FragmentCounter]);
    334161      } else {
    335162        DoeLog(1) && (eLog()<< Verbose(1) << "Subgraph " << MolecularWalker << " has no atoms!" << endl);
    … …  
    399226      // store Adjacency file
    400227      std::string filename = prefix + ADJACENCYFILE;
    401       StoreAdjacencyToFile(filename);
     228      mol->StoreAdjacencyToFile(filename);
    402229    }
    403230
    … …  
    434261
    435262
     263/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
     264 * -# constructs a complete keyset of the molecule
     265 * -# In a loop over all possible roots from the given rootstack
     266 *  -# increases order of root site
     267 *  -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
     268 *  -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
     269as the restricted one and each site in the set as the root)
     270 *  -# these are merged into a fragment list of keysets
     271 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
     272 * Important only is that we create all fragments, it is not important if we create them more than once
     273 * as these copies are filtered out via use of the hash table (KeySet).
     274 * \param *out output stream for debugging
     275 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
     276 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
     277 * \return pointer to Graph list
     278 */
     279void Fragmentation::FragmentBOSSANOVA(molecule *mol, Graph *&FragmentList, KeyStack &RootStack)
     280{
     281  Graph ***FragmentLowerOrdersList = NULL;
     282  int NumLevels = 0;
     283  int NumMolecules = 0;
     284  int TotalNumMolecules = 0;
     285  int *NumMoleculesOfOrder = NULL;
     286  int Order = 0;
     287  int UpgradeCount = RootStack.size();
     288  KeyStack FragmentRootStack;
     289  int RootKeyNr = 0;
     290  int RootNr = 0;
     291  UniqueFragments FragmentSearch;
     292
     293  DoLog(0) && (Log() << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl);
     294
     295  // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
     296  // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
     297  NumMoleculesOfOrder = new int[UpgradeCount];
     298  FragmentLowerOrdersList = new Graph**[UpgradeCount];
     299
     300  for(int i=0;i<UpgradeCount;i++) {
     301    NumMoleculesOfOrder[i] = 0;
     302    FragmentLowerOrdersList[i] = NULL;
     303  }
     304
     305  FragmentSearch.Init(mol->FindAtom(RootKeyNr));
     306
     307  // Construct the complete KeySet which we need for topmost level only (but for all Roots)
     308  KeySet CompleteMolecule;
     309  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     310    CompleteMolecule.insert((*iter)->GetTrueFather()->getNr());
     311  }
     312
     313  // 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
     314  // 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),
     315  // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
     316  // 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)
     317  RootNr = 0;   // counts through the roots in RootStack
     318  while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
     319    RootKeyNr = RootStack.front();
     320    RootStack.pop_front();
     321    atom *Walker = mol->FindAtom(RootKeyNr);
     322    // check cyclic lengths
     323    //if ((MinimumRingSize[Walker->GetTrueFather()->getNr()] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->getNr()])) {
     324    //  Log() << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
     325    //} else
     326    {
     327      // increase adaptive order by one
     328      Walker->GetTrueFather()->AdaptiveOrder++;
     329      Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
     330
     331      // initialise Order-dependent entries of UniqueFragments structure
     332      class PowerSetGenerator PSG(&FragmentSearch, Walker->AdaptiveOrder);
     333
     334      // allocate memory for all lower level orders in this 1D-array of ptrs
     335      NumLevels = 1 << (Order-1); // (int)pow(2,Order);
     336      FragmentLowerOrdersList[RootNr] = new Graph*[NumLevels];
     337      for (int i=0;i<NumLevels;i++)
     338        FragmentLowerOrdersList[RootNr][i] = NULL;
     339
     340      // create top order where nothing is reduced
     341      DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
     342      DoLog(0) && (Log() << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl); // , NumLevels is " << NumLevels << "
     343
     344      // Create list of Graphs of current Bond Order (i.e. F_{ij})
     345      FragmentLowerOrdersList[RootNr][0] = new Graph;
     346      FragmentSearch.PrepareForPowersetGeneration(1., FragmentLowerOrdersList[RootNr][0], Walker);
     347      NumMoleculesOfOrder[RootNr] = PSG(CompleteMolecule);
     348
     349      // output resulting number
     350      DoLog(1) && (Log() << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl);
     351      if (NumMoleculesOfOrder[RootNr] != 0) {
     352        NumMolecules = 0;
     353      } else {
     354        Walker->GetTrueFather()->MaxOrder = true;
     355      }
     356      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
     357      //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
     358      TotalNumMolecules += NumMoleculesOfOrder[RootNr];
     359//      Log() << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     360      RootStack.push_back(RootKeyNr); // put back on stack
     361      RootNr++;
     362    }
     363  }
     364  DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
     365  DoLog(1) && (Log() << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl);
     366  DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
     367
     368  // cleanup FragmentSearch structure
     369  FragmentSearch.Cleanup();
     370
     371  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
     372  // 5433222211111111
     373  // 43221111
     374  // 3211
     375  // 21
     376  // 1
     377
     378  // Subsequently, we combine all into a single list (FragmentList)
     379  CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, mol);
     380  FreeAllOrdersList(FragmentLowerOrdersList, RootStack, mol);
     381  delete[](NumMoleculesOfOrder);
     382
     383  DoLog(0) && (Log() << Verbose(0) << "End of FragmentBOSSANOVA." << endl);
     384};
     385
     386/** Stores a fragment from \a KeySet into \a molecule.
     387 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
     388 * molecule and adds missing hydrogen where bonds were cut.
     389 * \param *out output stream for debugging messages
     390 * \param &Leaflet pointer to KeySet structure
     391 * \param IsAngstroem whether we have Ansgtroem or bohrradius
     392 * \return pointer to constructed molecule
     393 */
     394molecule * Fragmentation::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
     395{
     396  atom **SonList = new atom*[mol->getAtomCount()];
     397  molecule *Leaf = World::getInstance().createMolecule();
     398
     399  for(int i=0;i<mol->getAtomCount();i++)
     400    SonList[i] = NULL;
     401
     402//  Log() << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
     403  StoreFragmentFromKeySet_Init(mol, Leaf, Leaflet, SonList);
     404  // create the bonds between all: Make it an induced subgraph and add hydrogen
     405//  Log() << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
     406  CreateInducedSubgraphOfFragment(mol, Leaf, SonList, IsAngstroem);
     407
     408  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
     409  delete[](SonList);
     410//  Log() << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
     411  return Leaf;
     412};
     413
     414
     415/** Estimates by educated guessing (using upper limit) the expected number of fragments.
     416 * The upper limit is
     417 * \f[
     418 *  n = N \cdot C^k
     419 * \f]
     420 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
     421 * \param *out output stream for debugging
     422 * \param order bond order k
     423 * \return number n of fragments
     424 */
     425int Fragmentation::GuesstimateFragmentCount(int order)
     426{
     427  size_t c = 0;
     428  int FragmentCount;
     429  // get maximum bond degree
     430  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     431    const BondList& ListOfBonds = (*iter)->getListOfBonds();
     432    c = (ListOfBonds.size() > c) ? ListOfBonds.size() : c;
     433  }
     434  FragmentCount = mol->NoNonHydrogen*(1 << (c*order));
     435  DoLog(1) && (Log() << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << mol->NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl);
     436  return FragmentCount;
     437};
     438
     439
     440/** Checks whether the OrderAtSite is still below \a Order at some site.
     441 * \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
     442 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
     443 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
     444 * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
     445 * \return true - needs further fragmentation, false - does not need fragmentation
     446 */
     447bool Fragmentation::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, std::string path)
     448{
     449  bool status = false;
     450
     451  // initialize mask list
     452  for(int i=mol->getAtomCount();i--;)
     453    AtomMask[i] = false;
     454
     455  if (Order < 0) { // adaptive increase of BondOrder per site
     456    if (AtomMask[mol->getAtomCount()] == true)  // break after one step
     457      return false;
     458
     459    // transmorph graph keyset list into indexed KeySetList
     460    if (GlobalKeySetList == NULL) {
     461      DoeLog(1) && (eLog()<< Verbose(1) << "Given global key set list (graph) is NULL!" << endl);
     462      return false;
     463    }
     464    map<int,KeySet> *IndexKeySetList = GraphToIndexedKeySet(GlobalKeySetList);
     465
     466    // parse the EnergyPerFragment file
     467    map<int, pair<double,int> > *AdaptiveCriteriaList = ScanAdaptiveFileIntoMap(path, *IndexKeySetList); // (Root No., (Value, Order)) !
     468    if (AdaptiveCriteriaList->empty()) {
     469      DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
     470      for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     471    #ifdef ADDHYDROGEN
     472        if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
     473    #endif
     474        {
     475          AtomMask[(*iter)->getNr()] = true;  // include all (non-hydrogen) atoms
     476          status = true;
     477        }
     478      }
     479    }
     480    // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
     481    map<double, pair<int,int> > *FinalRootCandidates = ReMapAdaptiveCriteriaListToValue(AdaptiveCriteriaList, mol);
     482
     483    // pick the ones still below threshold and mark as to be adaptively updated
     484    MarkUpdateCandidates(AtomMask, *FinalRootCandidates, Order, mol);
     485
     486    delete[](IndexKeySetList);
     487    delete[](AdaptiveCriteriaList);
     488    delete[](FinalRootCandidates);
     489  } else { // global increase of Bond Order
     490    for(molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     491  #ifdef ADDHYDROGEN
     492      if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
     493  #endif
     494      {
     495        AtomMask[(*iter)->getNr()] = true;  // include all (non-hydrogen) atoms
     496        if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->getNr()]))
     497          status = true;
     498      }
     499    }
     500    if ((!Order) && (!AtomMask[mol->getAtomCount()]))  // single stepping, just check
     501      status = true;
     502
     503    if (!status) {
     504      if (Order == 0)
     505        DoLog(1) && (Log() << Verbose(1) << "Single stepping done." << endl);
     506      else
     507        DoLog(1) && (Log() << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl);
     508    }
     509  }
     510
     511  PrintAtomMask(AtomMask, mol->getAtomCount()); // for debugging
     512
     513  return status;
     514};
     515
    436516/** Stores pairs (Atom::Nr, Atom::AdaptiveOrder) into file.
    437517 * Atoms not present in the file get "-1".
    … …  
    439519 * \return true - file writable, false - not writable
    440520 */
    441 bool molecule::StoreOrderAtSiteFile(std::string &path)
     521bool Fragmentation::StoreOrderAtSiteFile(std::string &path)
    442522{
    443523  string line;
    … …  
    448528  DoLog(1) && (Log() << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl);
    449529  if (file.good()) {
    450     for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputOrder), &file));
     530    for_each(mol->begin(),mol->end(),bind2nd(mem_fun(&atom::OutputOrder), &file));
    451531    file.close();
    452532    DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    … …  
    458538};
    459539
     540
    460541/** Parses pairs(Atom::Nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
    461542 * Atoms not present in the file get "0".
    … …  
    464545 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
    465546 */
    466 bool molecule::ParseOrderAtSiteFromFile(std::string &path)
    467 {
    468   unsigned char *OrderArray = new unsigned char[getAtomCount()];
    469   bool *MaxArray = new bool[getAtomCount()];
     547bool Fragmentation::ParseOrderAtSiteFromFile(std::string &path)
     548{
     549  unsigned char *OrderArray = new unsigned char[mol->getAtomCount()];
     550  bool *MaxArray = new bool[mol->getAtomCount()];
    470551  bool status;
    471552  int AtomNr, value;
    … …  
    473554  ifstream file;
    474555
    475   for(int i=0;i<getAtomCount();i++) {
     556  for(int i=0;i<mol->getAtomCount();i++) {
    476557    OrderArray[i] = 0;
    477558    MaxArray[i] = false;
    … …  
    496577
    497578    // set atom values
    498     for(internal_iterator iter=atoms.begin();iter!=atoms.end();++iter){
     579    for(molecule::const_iterator iter=mol->begin();iter!=mol->end();++iter){
    499580      (*iter)->AdaptiveOrder = OrderArray[(*iter)->getNr()];
    500581      (*iter)->MaxOrder = MaxArray[(*iter)->getNr()];
    … …  
    514595  DoLog(1) && (Log() << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl);
    515596  return status;
     597};
     598
     599/** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
     600 * \param *out output stream for debugging
     601 * \param *&SortIndex Mapping array of size molecule::AtomCount
     602 * \return true - success, false - failure of SortIndex alloc
     603 */
     604bool Fragmentation::CreateMappingLabelsToConfigSequence(int *&SortIndex)
     605{
     606  if (SortIndex != NULL) {
     607    DoLog(1) && (Log() << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl);
     608    return false;
     609  }
     610  SortIndex = new int[mol->getAtomCount()];
     611  for(int i=mol->getAtomCount();i--;)
     612    SortIndex[i] = -1;
     613
     614  int AtomNo = 0;
     615  for(molecule::const_iterator iter=mol->begin();iter!=mol->end();++iter){
     616    ASSERT(SortIndex[(*iter)->getNr()]==-1,"Same SortIndex set twice");
     617    SortIndex[(*iter)->getNr()] = AtomNo++;
     618  }
     619
     620  return true;
    516621};
    517622
    … …  
    524629 * \return number of atoms in fragment
    525630 */
    526 int StoreFragmentFromKeySet_Init(molecule *mol, molecule *Leaf, KeySet &Leaflet, atom **SonList)
     631int Fragmentation::StoreFragmentFromKeySet_Init(molecule *mol, molecule *Leaf, KeySet &Leaflet, atom **SonList)
    527632{
    528633  atom *FatherOfRunner = NULL;
    … …  
    538643};
    539644
     645
    540646/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
    541647 * \param *out output stream for debugging messages
    … …  
    545651 * \param **SonList list which atom of \a *Leaf is a son of which atom in \a *mol
    546652 */
    547 void CreateInducedSubgraphOfFragment(molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
     653void Fragmentation::CreateInducedSubgraphOfFragment(molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
    548654{
    549655  bool LonelyFlag = false;
    … …  
    604710};
    605711
    606 /** Stores a fragment from \a KeySet into \a molecule.
    607  * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
    608  * molecule and adds missing hydrogen where bonds were cut.
    609  * \param *out output stream for debugging messages
    610  * \param &Leaflet pointer to KeySet structure
    611  * \param IsAngstroem whether we have Ansgtroem or bohrradius
    612  * \return pointer to constructed molecule
    613  */
    614 molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
    615 {
    616   atom **SonList = new atom*[getAtomCount()];
    617   molecule *Leaf = World::getInstance().createMolecule();
    618 
    619   for(int i=0;i<getAtomCount();i++)
    620     SonList[i] = NULL;
    621 
    622 //  Log() << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
    623   StoreFragmentFromKeySet_Init(this, Leaf, Leaflet, SonList);
    624   // create the bonds between all: Make it an induced subgraph and add hydrogen
    625 //  Log() << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
    626   CreateInducedSubgraphOfFragment(this, Leaf, SonList, IsAngstroem);
    627 
    628   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
    629   delete[](SonList);
    630 //  Log() << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
    631   return Leaf;
    632 };
    633 
    634 /** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
    635  * -# constructs a complete keyset of the molecule
    636  * -# In a loop over all possible roots from the given rootstack
    637  *  -# increases order of root site
    638  *  -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
    639  *  -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
    640 as the restricted one and each site in the set as the root)
    641  *  -# these are merged into a fragment list of keysets
    642  * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
    643  * Important only is that we create all fragments, it is not important if we create them more than once
    644  * as these copies are filtered out via use of the hash table (KeySet).
    645  * \param *out output stream for debugging
    646  * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
    647  * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
    648  * \return pointer to Graph list
    649  */
    650 void molecule::FragmentBOSSANOVA(Graph *&FragmentList, KeyStack &RootStack)
    651 {
    652   Graph ***FragmentLowerOrdersList = NULL;
    653   int NumLevels = 0;
    654   int NumMolecules = 0;
    655   int TotalNumMolecules = 0;
    656   int *NumMoleculesOfOrder = NULL;
    657   int Order = 0;
    658   int UpgradeCount = RootStack.size();
    659   KeyStack FragmentRootStack;
    660   int RootKeyNr = 0;
    661   int RootNr = 0;
    662   struct UniqueFragments FragmentSearch;
    663 
    664   DoLog(0) && (Log() << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl);
    665 
    666   // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
    667   // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
    668   NumMoleculesOfOrder = new int[UpgradeCount];
    669   FragmentLowerOrdersList = new Graph**[UpgradeCount];
    670 
    671   for(int i=0;i<UpgradeCount;i++) {
    672     NumMoleculesOfOrder[i] = 0;
    673     FragmentLowerOrdersList[i] = NULL;
    674   }
    675 
    676   FragmentSearch.Init(FindAtom(RootKeyNr));
    677 
    678   // Construct the complete KeySet which we need for topmost level only (but for all Roots)
    679   KeySet CompleteMolecule;
    680   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    681     CompleteMolecule.insert((*iter)->GetTrueFather()->getNr());
    682   }
    683 
    684   // 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
    685   // 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),
    686   // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
    687   // 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)
    688   RootNr = 0;   // counts through the roots in RootStack
    689   while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
    690     RootKeyNr = RootStack.front();
    691     RootStack.pop_front();
    692     atom *Walker = FindAtom(RootKeyNr);
    693     // check cyclic lengths
    694     //if ((MinimumRingSize[Walker->GetTrueFather()->getNr()] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->getNr()])) {
    695     //  Log() << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
    696     //} else
    697     {
    698       // increase adaptive order by one
    699       Walker->GetTrueFather()->AdaptiveOrder++;
    700       Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
    701 
    702       // initialise Order-dependent entries of UniqueFragments structure
    703       class PowerSetGenerator PSG(&FragmentSearch, Walker->AdaptiveOrder);
    704 
    705       // allocate memory for all lower level orders in this 1D-array of ptrs
    706       NumLevels = 1 << (Order-1); // (int)pow(2,Order);
    707       FragmentLowerOrdersList[RootNr] = new Graph*[NumLevels];
    708       for (int i=0;i<NumLevels;i++)
    709         FragmentLowerOrdersList[RootNr][i] = NULL;
    710 
    711       // create top order where nothing is reduced
    712       DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
    713       DoLog(0) && (Log() << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl); // , NumLevels is " << NumLevels << "
    714 
    715       // Create list of Graphs of current Bond Order (i.e. F_{ij})
    716       FragmentLowerOrdersList[RootNr][0] = new Graph;
    717       FragmentSearch.PrepareForPowersetGeneration(1., FragmentLowerOrdersList[RootNr][0], Walker);
    718       NumMoleculesOfOrder[RootNr] = PSG(CompleteMolecule);
    719 
    720       // output resulting number
    721       DoLog(1) && (Log() << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl);
    722       if (NumMoleculesOfOrder[RootNr] != 0) {
    723         NumMolecules = 0;
    724       } else {
    725         Walker->GetTrueFather()->MaxOrder = true;
    726       }
    727       // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
    728       //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
    729       TotalNumMolecules += NumMoleculesOfOrder[RootNr];
    730 //      Log() << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    731       RootStack.push_back(RootKeyNr); // put back on stack
    732       RootNr++;
    733     }
    734   }
    735   DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
    736   DoLog(1) && (Log() << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl);
    737   DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
    738 
    739   // cleanup FragmentSearch structure
    740   FragmentSearch.Cleanup();
    741 
    742   // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
    743   // 5433222211111111
    744   // 43221111
    745   // 3211
    746   // 21
    747   // 1
    748 
    749   // Subsequently, we combine all into a single list (FragmentList)
    750   CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, this);
    751   FreeAllOrdersList(FragmentLowerOrdersList, RootStack, this);
    752   delete[](NumMoleculesOfOrder);
    753 
    754   DoLog(0) && (Log() << Verbose(0) << "End of FragmentBOSSANOVA." << endl);
    755 };
    756 
    757 /** Corrects the nuclei position if the fragment was created over the cell borders.
    758  * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
    759  * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
    760  * and re-add the bond. Looping on the distance check.
    761  * \param *out ofstream for debugging messages
    762  */
    763 bool molecule::ScanForPeriodicCorrection()
    764 {
    765   bond *Binder = NULL;
    766   //bond *OtherBinder = NULL;
    767   atom *Walker = NULL;
    768   atom *OtherWalker = NULL;
    769   RealSpaceMatrix matrix = World::getInstance().getDomain().getM();
    770   enum GraphEdge::Shading *ColorList = NULL;
    771   double tmp;
    772   //bool LastBond = true; // only needed to due list construct
    773   Vector Translationvector;
    774   //std::deque<atom *> *CompStack = NULL;
    775   std::deque<atom *> *AtomStack = new std::deque<atom *>; // (getAtomCount());
    776   bool flag = true;
    777   BondGraph *BG = World::getInstance().getBondGraph();
    778 
    779   DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl);
    780 
    781   ColorList = new enum GraphEdge::Shading[getAtomCount()];
    782   for (int i=0;i<getAtomCount();i++)
    783     ColorList[i] = (enum GraphEdge::Shading)0;
    784   if (flag) {
    785     // remove bonds that are beyond bonddistance
    786     Translationvector.Zero();
    787     // scan all bonds
    788     flag = false;
    789     for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner) {
    790       const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    791       for(BondList::const_iterator BondRunner = ListOfBonds.begin();
    792           (!flag) && (BondRunner != ListOfBonds.end());
    793           ++BondRunner) {
    794         Binder = (*BondRunner);
    795         for (int i=NDIM;i--;) {
    796           tmp = fabs(Binder->leftatom->at(i) - Binder->rightatom->at(i));
    797           //Log() << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
    798           const range<double> MinMaxDistance(
    799               BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
    800           if (!MinMaxDistance.isInRange(tmp)) {
    801             DoLog(2) && (Log() << Verbose(2) << "Correcting at bond " << *Binder << "." << endl);
    802             flag = true;
    803             break;
    804           }
    805         }
    806       }
    807     }
    808     //if (flag) {
    809     if (0) {
    810       // create translation vector from their periodically modified distance
    811       for (int i=NDIM;i--;) {
    812         tmp = Binder->leftatom->at(i) - Binder->rightatom->at(i);
    813         const range<double> MinMaxDistance(
    814             BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
    815         if (fabs(tmp) > MinMaxDistance.last)  // check against Min is not useful for components
    816           Translationvector[i] = (tmp < 0) ? +1. : -1.;
    817       }
    818       Translationvector *= matrix;
    819       //Log() << Verbose(3) << "Translation vector is ";
    820       Log() << Verbose(0) << Translationvector <<  endl;
    821       // apply to all atoms of first component via BFS
    822       for (int i=getAtomCount();i--;)
    823         ColorList[i] = GraphEdge::white;
    824       AtomStack->push_front(Binder->leftatom);
    825       while (!AtomStack->empty()) {
    826         Walker = AtomStack->front();
    827         AtomStack->pop_front();
    828         //Log() << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
    829         ColorList[Walker->getNr()] = GraphEdge::black;    // mark as explored
    830         *Walker += Translationvector; // translate
    831         const BondList& ListOfBonds = Walker->getListOfBonds();
    832         for (BondList::const_iterator Runner = ListOfBonds.begin();
    833             Runner != ListOfBonds.end();
    834             ++Runner) {
    835           if ((*Runner) != Binder) {
    836             OtherWalker = (*Runner)->GetOtherAtom(Walker);
    837             if (ColorList[OtherWalker->getNr()] == GraphEdge::white) {
    838               AtomStack->push_front(OtherWalker); // push if yet unexplored
    839             }
    840           }
    841         }
    842       }
    843 //      // re-add bond
    844 //      if (OtherBinder == NULL) { // is the only bond?
    845 //        //Do nothing
    846 //      } else {
    847 //        if (!LastBond) {
    848 //          link(Binder, OtherBinder); // no more implemented bond::previous ...
    849 //        } else {
    850 //          link(OtherBinder, Binder); // no more implemented bond::previous ...
    851 //        }
    852 //      }
    853     } else {
    854       DoLog(3) && (Log() << Verbose(3) << "No corrections for this fragment." << endl);
    855     }
    856     //delete(CompStack);
    857   }
    858   // free allocated space from ReturnFullMatrixforSymmetric()
    859   delete(AtomStack);
    860   delete[](ColorList);
    861   DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
    862 
    863   return flag;
    864 };

  • src/Fragmentation/Makefile.am

    r212c179 r246e13  
    66        Fragmentation/EnergyMatrix.cpp \
    77        Fragmentation/ForceMatrix.cpp \
     8        Fragmentation/Fragmentation.cpp \
    89        Fragmentation/fragmentation_helpers.cpp \
    910        Fragmentation/helpers.cpp \
    … …  
    1920        Fragmentation/EnergyMatrix.hpp \
    2021        Fragmentation/ForceMatrix.hpp \
     22        Fragmentation/Fragmentation.hpp \
    2123        Fragmentation/fragmentation_helpers.hpp \
    2224        Fragmentation/helpers.hpp \
Note: See TracChangeset for help on using the changeset viewer.