Changeset 4aa03a for src/molecules.cpp


Ignore:
Timestamp:
Jul 31, 2008, 2:05:37 PM (16 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:
fb9364
Parents:
af6d8f
Message:

new test function CheckForConnectedSubgraph() and bigger rewrite of PowerSetGenerator() and SPFragmentGenerator() to allow for ings

CheckForConnectedSubgraph(): Is an O(N2) function to check whether we do not create fragments we don't want to have (i.e. disconnected ones)
PowerSetGenerator(): Instead of adding the possible Walkers to an AtomStack, we use the bonds we store in the BondsPerSPList, which by the way have the necessary directional information (i.e. they come with the predecessor in Bond#leftatom). This way, we generate fragments in ring structures without any problems.
SPFragmentGenerator(): had to be adapted as the calling in PowerSetGenerator() had changed a bit. We already let the root be added by this function and don't add it in PowerSetGenerator(), which simplifies (no if-catch for first order) the code there a bit. Furthermore, we had to check whether an added vertex was not already in the keyset (this may happen, as both edge 1<->2 as 2<->1 may appear at different SPLevels in BondsPerSPLevel: Ring structures in contrast to tree structures don not allow for uniqueness! This is first realised through the hash map storing the keysets). This again caused a change in recognizing when the stack contains a fragment, as Suborder-bits (bits is the number of atoms added for this combination) is not the way anymore - we use an additional Counter "Added".

Test on Benzonitrile was successful. We counted the possible fragments by hand, checked them furthermore till order 5 and the numbers were ok up till full order 8 in the end.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecules.cpp

    raf6d8f r4aa03a  
    25232523        No = (*runner).second.first;
    25242524        Walker = FindAtom(No);
    2525         if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
     2525        //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
    25262526          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
    25272527          AtomMask[No] = true;
    25282528          status = true;
    2529         } else
    2530           *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
     2529        //} else
     2530          //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
    25312531      }
    25322532      // close and done
    … …  
    25562556      {
    25572557        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
    2558         if ((Order != 0) && (Walker->AdaptiveOrder < Order) && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
     2558        if ((Order != 0) && (Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
    25592559          status = true;
    25602560      }
    25612561    }
    25622562    if ((Order == 0) && (AtomMask[AtomCount] == true))  // single stepping, just check
    2563       status = false;
     2563      status = true;
    25642564     
    25652565    if (!status)
    … …  
    31523152  atom **SonList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::StoreFragmentFromStack: **SonList");
    31533153  molecule *Leaf = new molecule(elemente);
     3154  bool LonelyFlag = false;
     3155  int size;
    31543156 
    31553157//  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
    … …  
    31643166
    31653167  // first create the minimal set of atoms from the KeySet
     3168  size = 0;
    31663169  for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
    31673170    FatherOfRunner = FindAtom((*runner));  // find the id
    31683171    SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
     3172    size++;
    31693173  }
    31703174 
    … …  
    31743178  while (Runner->next != Leaf->end) {
    31753179    Runner = Runner->next;
     3180    LonelyFlag = true;
    31763181    FatherOfRunner = Runner->father;
    31773182    if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
    … …  
    31913196//            *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
    31923197          }
     3198          LonelyFlag = false;
    31933199        } else {
    31943200//          *out << ", who has no son in this fragment molecule." << endl;
    … …  
    32023208    } else {
    32033209      *out << Verbose(0) << "ERROR: Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
     3210    }
     3211    if ((LonelyFlag) && (size > 1)) {
     3212      *out << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl;
    32043213    }
    32053214#ifdef ADDHYDROGEN
    … …  
    34683477  int NumCombinations;
    34693478  bool bit;
    3470   int bits, TouchedIndex, SubSetDimension, SP;
     3479  int bits, TouchedIndex, SubSetDimension, SP, Added;
    34713480  int Removal;
     3481  int SpaceLeft;
    34723482  int *TouchedList = (int *) Malloc(sizeof(int)*(SubOrder+1), "molecule::SPFragmentGenerator: *TouchedList");
    34733483  bond *Binder = NULL;
    34743484  bond **BondsList = NULL;
     3485  KeySetTestPair TestKeySetInsert;
    34753486
    34763487  NumCombinations = 1 << SetDimension;
    … …  
    35003511    if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
    35013512      // --1-- add this set of the power set of bond partners to the snake stack
     3513      Added = 0;
    35023514      for (int j=0;j<SetDimension;j++) {  // pull out every bit by shifting
    35033515        bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
    … …  
    35083520            //*out << Verbose(2+verbosity) << "Not used so far, label " << FragmentSearch->Labels[OtherWalker->nr] << " is bigger than Root's " << FragmentSearch->Labels[FragmentSearch->Root->nr] << "." << endl;
    35093521          *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
    3510           TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
    3511           FragmentSearch->FragmentSet->insert(OtherWalker->nr);
     3522          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
     3523          if (TestKeySetInsert.second) {
     3524            TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
     3525            Added++;
     3526          } else {
     3527            *out << Verbose(2+verbosity) << "This was item was already present in the keyset." << endl;
     3528          }
    35123529            //FragmentSearch->UsedList[OtherWalker->nr][i] = true;
    35133530          //}
    … …  
    35173534      }
    35183535     
    3519       if (bits < SubOrder) {
    3520         *out << Verbose(1+verbosity) << "There's still some space left on stack: " << (SubOrder - bits) << "." << endl;
    3521         // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
    3522         SP = RootDistance+1;    // this is the next level
    3523         // first count the members in the subset
    3524         SubSetDimension = 0;
    3525         Binder = FragmentSearch->BondsPerSPList[2*SP];          // start node for this level
    3526         while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {                // compare to end node of this level
    3527           Binder = Binder->next;
    3528           for (int k=TouchedIndex;k--;) {
    3529             if (Binder->Contains(TouchedList[k]))   // if we added this very endpiece
    3530               SubSetDimension++;
     3536      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
     3537      if (SpaceLeft > 0) {
     3538        *out << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
     3539        if (SubOrder > 1) {    // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
     3540          // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
     3541          SP = RootDistance+1;  // this is the next level
     3542          // first count the members in the subset
     3543          SubSetDimension = 0;
     3544          Binder = FragmentSearch->BondsPerSPList[2*SP];                // start node for this level
     3545          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {              // compare to end node of this level
     3546            Binder = Binder->next;
     3547            for (int k=TouchedIndex;k--;) {
     3548              if (Binder->Contains(TouchedList[k]))   // if we added this very endpiece
     3549                SubSetDimension++;
     3550            }
    35313551          }
     3552          // then allocate and fill the list
     3553          BondsList = (bond **) Malloc(sizeof(bond *)*SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
     3554          SubSetDimension = 0;
     3555          Binder = FragmentSearch->BondsPerSPList[2*SP];
     3556          while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {
     3557            Binder = Binder->next;
     3558            for (int k=0;k<TouchedIndex;k++) {
     3559              if (Binder->leftatom->nr == TouchedList[k])   // leftatom is always the close one
     3560                BondsList[SubSetDimension++] = Binder;
     3561            }
     3562          }
     3563          *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
     3564          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits); 
     3565          Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
    35323566        }
    3533         // then allocate and fill the list
    3534         BondsList = (bond **) Malloc(sizeof(bond *)*SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
    3535         SubSetDimension = 0;
    3536         Binder = FragmentSearch->BondsPerSPList[2*SP];
    3537         while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {
    3538           Binder = Binder->next;
    3539           for (int k=0;k<TouchedIndex;k++) {
    3540             if (Binder->leftatom->nr == TouchedList[k])   // leftatom is always the close one
    3541               BondsList[SubSetDimension++] = Binder;
    3542           }
    3543         }
    3544         *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
    3545         SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits); 
    3546         Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
    35473567      } else {
    35483568        // --2-- otherwise store the complete fragment
    … …  
    35523572        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
    35533573          *out << (*runner) << " ";
     3574        *out << endl;
     3575        if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
     3576          *out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl;
    35543577        InsertFragmentIntoGraph(out, FragmentSearch);
    35553578        //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
    … …  
    35773600  *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
    35783601};
     3602
     3603/** For a given keyset \a *Fragment, checks whether it is connected in the current molecule.
     3604 * \param *out output stream for debugging
     3605 * \param *Fragment Keyset of fragment's vertices
     3606 * \return true - connected, false - disconnected
     3607 * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
     3608 */
     3609bool molecule::CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment)
     3610{
     3611  atom *Walker = NULL, *Walker2 = NULL;
     3612  bool BondStatus = false;
     3613  int size;
     3614 
     3615  *out << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl;
     3616  *out << Verbose(2) << "Disconnected atom: ";
     3617 
     3618  // count number of atoms in graph
     3619  size = 0;
     3620  for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
     3621    size++;
     3622  if (size > 1)
     3623    for(KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++) {
     3624      Walker = FindAtom(*runner);
     3625      BondStatus = false;
     3626      for(KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
     3627        Walker2 = FindAtom(*runners);
     3628        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
     3629          if (ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker) == Walker2) {
     3630            BondStatus = true;
     3631            break;
     3632          }
     3633          if (BondStatus)
     3634            break;
     3635        }
     3636      }
     3637      if (!BondStatus) {
     3638        *out << (*Walker) << endl;
     3639        return false;
     3640      }
     3641    }
     3642  else {
     3643    *out << "none." << endl;
     3644    return true;
     3645  }
     3646  *out << "none." << endl;
     3647
     3648  *out << Verbose(1) << "End of CheckForConnectedSubgraph" << endl;
     3649
     3650  return true;
     3651}
    35793652
    35803653/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
    … …  
    35963669{
    35973670  int SP, UniqueIndex, AtomKeyNr;
    3598   int *NumberOfAtomsSPLevel = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *SPLevelCount");
    3599   atom *Walker = NULL, *OtherWalker = NULL;
     3671  atom *Walker = NULL, *OtherWalker = NULL, *Predecessor = NULL;
    36003672  bond *Binder = NULL;
     3673  bond *CurrentEdge = NULL;
    36013674  bond **BondsList = NULL;
    3602   KeyStack AtomStack;
    3603   atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
    36043675  int RootKeyNr = FragmentSearch.Root->nr;
    36053676  int Counter = FragmentSearch.FragmentCounter;
    3606 
    3607   for (int i=AtomCount;i--;) {
    3608     PredecessorList[i] = NULL;
    3609   }
     3677  int RemainingWalkers;
    36103678
    36113679  *out << endl;
    36123680  *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;
    36133681
     3682  // prepare Label and SP arrays of the BFS search
    36143683  UniqueIndex = 0;
    36153684  if (FragmentSearch.Labels[RootKeyNr] == -1)
    36163685    FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
    36173686  FragmentSearch.ShortestPathList[RootKeyNr] = 0;
    3618   // prepare the atom stack counters (number of atoms with certain SP on stack)
    3619   for (int i=Order;i--;)
    3620     NumberOfAtomsSPLevel[i] = 0;
    3621   NumberOfAtomsSPLevel[0] = 1;  // for root
    3622   SP = -1;
     3687
     3688  // prepare root level (SP = 0) and a loop bond denoting Root
     3689  for (int i=1;i<Order;i++)
     3690    FragmentSearch.BondsPerSPCount[i] = 0;
     3691  FragmentSearch.BondsPerSPCount[0] = 1;
     3692  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
     3693  add(Binder, FragmentSearch.BondsPerSPList[1]);
     3694 
     3695  // do a BFS search to fill the SP lists and label the found vertices
     3696  // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
     3697  // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
     3698  // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
     3699  // (EdgeinSPLevel) of this tree ...
     3700  // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
     3701  // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
    36233702  *out << endl;
    36243703  *out << Verbose(0) << "Starting BFS analysis ..." << endl;
    3625   // push as first on atom stack and goooo ...
    3626   AtomStack.clear(); 
    3627   AtomStack.push_back(RootKeyNr);
    3628   *out << Verbose(0) << "Cleared AtomStack and pushed root as first item onto it." << endl;
    3629   // do a BFS search to fill the SP lists and label the found vertices
    3630   while (!AtomStack.empty()) {
    3631     // pop next atom
    3632     AtomKeyNr = AtomStack.front();
    3633     AtomStack.pop_front();
    3634     if (SP != -1)
    3635       NumberOfAtomsSPLevel[SP]--;
    3636     if ((SP == -1) || (NumberOfAtomsSPLevel[SP] == -1)) {
    3637       SP++;
    3638       NumberOfAtomsSPLevel[SP]--; // carry over "-1" to next level
    3639       *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with 0 item(s)";
    3640       if (SP > 0)
    3641         *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
    3642       else
    3643         *out << "." << endl;
     3704  for (SP = 0; SP < (Order-1); SP++) {
     3705    *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with " << FragmentSearch.BondsPerSPCount[SP] << " item(s)";
     3706    if (SP > 0) {
     3707      *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
    36443708      FragmentSearch.BondsPerSPCount[SP] = 0;
    3645     } else {
    3646       *out << Verbose(1) << "Still " << NumberOfAtomsSPLevel[SP]+1 << " on this SP (" << SP << ") level, currently having " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
    3647     }
    3648     Walker = FindAtom(AtomKeyNr);
    3649     *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << "." << endl;
    3650     // check for new sp level
    3651     // go through all its bonds
    3652     *out << Verbose(1) << "Going through all bonds of Walker." << endl;
    3653     for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
    3654       Binder = ListOfBondsPerAtom[AtomKeyNr][i];
    3655       OtherWalker = Binder->GetOtherAtom(Walker);
    3656       if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
    3657 #ifdef ADDHYDROGEN
    3658        && (OtherWalker->type->Z != 1)
    3659 #endif
    3660                                                             ) {  // skip hydrogens and restrict to fragment
    3661         *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
    3662         // set the label if not set (and push on root stack as well)
    3663         if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
    3664           FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
    3665           *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
    3666         } else {
    3667           *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
    3668         }
    3669         if ((OtherWalker != PredecessorList[AtomKeyNr]) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
    3670           FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
    3671           *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
    3672           if (FragmentSearch.ShortestPathList[OtherWalker->nr] < Order) { // only pass through those within reach of Order
    3673             // push for exploration on stack (only if the SP of OtherWalker is longer than of Walker! (otherwise it has been added already!)
    3674             if (FragmentSearch.ShortestPathList[OtherWalker->nr] > SP) {
    3675               if (FragmentSearch.ShortestPathList[OtherWalker->nr] < (Order-1)) {
    3676                 *out << Verbose(3) << "Putting on atom stack for further exploration." << endl;
    3677                 PredecessorList[OtherWalker->nr] = Walker;    // note down the one further up the exploration tree
    3678                 AtomStack.push_back(OtherWalker->nr);
    3679                 NumberOfAtomsSPLevel[FragmentSearch.ShortestPathList[OtherWalker->nr]]++;
    3680               } else {
    3681                 *out << Verbose(3) << "Not putting on atom stack, is at end of reach." << endl;
    3682               }
    3683               // add the bond in between to the SP list
    3684               Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
    3685               add(Binder, FragmentSearch.BondsPerSPList[2*SP+1]);
    3686               FragmentSearch.BondsPerSPCount[SP]++;
    3687               *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
    3688             } else *out << Verbose(3) << "Not putting on atom stack, nor adding bond, as " << *OtherWalker << "'s SP is shorter than Walker's." << endl;
    3689           } else *out << Verbose(3) << "Not continuing, as " << *OtherWalker << " is out of reach of order " << Order << "." << endl;
    3690         } else *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->nr << " is smaller than that of Root " << RootKeyNr << " or this is my predecessor." << endl;
    3691       } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
    3692     }
    3693   }
    3694   // reset predecessor list
    3695   for(int i=0;i<Order;i++) {
    3696     Binder = FragmentSearch.BondsPerSPList[2*i];
    3697     *out << Verbose(1) << "Current SP level is " << i << "." << endl;
    3698     while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
    3699       Binder = Binder->next;
    3700       PredecessorList[Binder->rightatom->nr] = NULL;  // By construction "OtherAtom" is always Bond::rightatom!
    3701     }
    3702   }
    3703   *out << endl;
     3709    } else
     3710      *out << "." << endl;
     3711
     3712    RemainingWalkers = FragmentSearch.BondsPerSPCount[SP];
     3713    CurrentEdge = FragmentSearch.BondsPerSPList[2*SP];    /// start of this SP level's list
     3714    while (CurrentEdge->next != FragmentSearch.BondsPerSPList[2*SP+1]) {    /// end of this SP level's list
     3715      CurrentEdge = CurrentEdge->next;
     3716      RemainingWalkers--;
     3717      Walker = CurrentEdge->rightatom;    // rightatom is always the one more distant
     3718      Predecessor = CurrentEdge->leftatom;    // ... and leftatom is predecessor
     3719      AtomKeyNr = Walker->nr;
     3720      *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
     3721      // check for new sp level
     3722      // go through all its bonds
     3723      *out << Verbose(1) << "Going through all bonds of Walker." << endl;
     3724      for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
     3725        Binder = ListOfBondsPerAtom[AtomKeyNr][i];
     3726        OtherWalker = Binder->GetOtherAtom(Walker);
     3727        if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
     3728  #ifdef ADDHYDROGEN
     3729         && (OtherWalker->type->Z != 1)
     3730  #endif
     3731                                                              ) {  // skip hydrogens and restrict to fragment
     3732          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
     3733          // set the label if not set (and push on root stack as well)
     3734          if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
     3735            FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
     3736            *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
     3737          } else {
     3738            *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
     3739          }
     3740          if ((OtherWalker != Predecessor) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
     3741            FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
     3742            *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
     3743            // add the bond in between to the SP list
     3744            Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
     3745            add(Binder, FragmentSearch.BondsPerSPList[2*(SP+1)+1]);
     3746            FragmentSearch.BondsPerSPCount[SP+1]++;
     3747            *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl;
     3748          } else *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->nr << " is smaller than that of Root " << RootKeyNr << " or this is my predecessor " << *Predecessor << "." << endl;
     3749        } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
     3750      }
     3751    }
     3752  }
    37043753 
    37053754  // outputting all list for debugging
    37063755  *out << Verbose(0) << "Printing all found lists." << endl;
    3707   for(int i=0;i<Order;i++) {
     3756  for(int i=1;i<Order;i++) {    // skip the root edge in the printing
    37083757    Binder = FragmentSearch.BondsPerSPList[2*i];
    37093758    *out << Verbose(1) << "Current SP level is " << i << "." << endl;
    … …  
    37153764 
    37163765  // creating fragments with the found edge sets  (may be done in reverse order, faster)
    3717   SP = 0;
     3766  SP = -1;  // the Root <-> Root edge must be subtracted!
    37183767  for(int i=Order;i--;) { // sum up all found edges
    37193768    Binder = FragmentSearch.BondsPerSPList[2*i];
    … …  
    37283777    *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
    37293778    FragmentSearch.FragmentSet->clear(); 
    3730     FragmentSearch.FragmentSet->insert(FragmentSearch.Root->nr);
     3779    *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
     3780    // prepare the subset and call the generator
     3781    BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
     3782    BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond
    37313783   
    3732     if (FragmentSearch.FragmentSet->size() == (unsigned int) Order) {
    3733       *out << Verbose(0) << "Enough items on stack already for a fragment!" << endl;
    3734       // store fragment as a KeySet
    3735       *out << Verbose(2) << "Found a new fragment[" << FragmentSearch.FragmentCounter << "], local nr.s are: ";
    3736       for(KeySet::iterator runner = FragmentSearch.FragmentSet->begin(); runner != FragmentSearch.FragmentSet->end(); runner++) {
    3737         *out << (*runner) << " ";
    3738       }
    3739       *out << endl;
    3740       InsertFragmentIntoGraph(out, &FragmentSearch);
    3741     } else {
    3742       *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
    3743       // prepare the subset and call the generator
    3744       BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
    3745       Binder = FragmentSearch.BondsPerSPList[0];
    3746       for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++) {
    3747         Binder = Binder->next;
    3748         BondsList[i] = Binder;
    3749       }
    3750       SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order-1); 
    3751       Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
    3752     }
     3784    SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
     3785   
     3786    Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
    37533787  } else {
    37543788    *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
    … …  
    37773811  }
    37783812
    3779   // free allocated memory
    3780   Free((void **)&NumberOfAtomsSPLevel, "molecule::PowerSetGenerator: *SPLevelCount");
    3781   Free((void **)&PredecessorList, "molecule::PowerSetGenerator: **PredecessorList");
    3782  
    37833813  // return list 
    37843814  *out << Verbose(0) << "End of PowerSetGenerator." << endl;
    … …  
    40394069    Walker = FindAtom(RootKeyNr);
    40404070    // check cyclic lengths
    4041     if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
    4042       *out << 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;
    4043     } else {
     4071    //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
     4072    //  *out << 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;
     4073    //} else
     4074    {
    40444075      // increase adaptive order by one
    40454076      Walker->GetTrueFather()->AdaptiveOrder++;
    … …  
    40604091      NumLevels = 1 << (Order-1); // (int)pow(2,Order);
    40614092      FragmentLowerOrdersList[RootNr] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
     4093      for (int i=0;i<NumLevels;i++)
     4094        FragmentLowerOrdersList[RootNr][i] = NULL;
    40624095     
    40634096      // create top order where nothing is reduced
    40644097      *out << Verbose(0) << "==============================================================================================================" << endl;
    4065       *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr-1) << " Roots remaining." << endl;
     4098      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl;
    40664099 
    40674100      // Create list of Graphs of current Bond Order (i.e. F_{ij})
    … …  
    40724105      NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
    40734106      *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    4074       NumMolecules = 0;
    4075      
    4076       if ((NumLevels >> 1) > 0) {
    4077         // create lower order fragments
    4078         *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
    4079         Order = Walker->AdaptiveOrder;
    4080         for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
    4081           // step down to next order at (virtual) boundary of powers of 2 in array
    4082           while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
    4083             Order--;
    4084           *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
    4085           for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
    4086             int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
    4087             *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
    4088             *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
    4089      
    4090             // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
    4091             //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
    4092             //NumMolecules = 0;
    4093             FragmentLowerOrdersList[RootNr][dest] = new Graph;
    4094             for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
    4095               for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    4096                 Graph TempFragmentList;
    4097                 FragmentSearch.TEFactor = -(*runner).second.second;
    4098                 FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
    4099                 FragmentSearch.Root = FindAtom(*sprinter);
    4100                 NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
    4101                 // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
    4102                 *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
    4103                 InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
     4107      if (NumMoleculesOfOrder[RootNr] != 0) {
     4108        NumMolecules = 0;
     4109       
     4110        if ((NumLevels >> 1) > 0) {
     4111          // create lower order fragments
     4112          *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
     4113          Order = Walker->AdaptiveOrder;
     4114          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
     4115            // step down to next order at (virtual) boundary of powers of 2 in array
     4116            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
     4117              Order--;
     4118            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
     4119            for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
     4120              int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
     4121              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
     4122              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
     4123       
     4124              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
     4125              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
     4126              //NumMolecules = 0;
     4127              FragmentLowerOrdersList[RootNr][dest] = new Graph;
     4128              for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
     4129                for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
     4130                  Graph TempFragmentList;
     4131                  FragmentSearch.TEFactor = -(*runner).second.second;
     4132                  FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
     4133                  FragmentSearch.Root = FindAtom(*sprinter);
     4134                  NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
     4135                  // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
     4136                  *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
     4137                  InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
     4138                }
    41044139              }
     4140              *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
    41054141            }
    4106             *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
    41074142          }
    41084143        }
     4144      } else {
     4145        *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
    41094146      }
    41104147      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
Note: See TracChangeset for help on using the changeset viewer.