Changeset fc850d for src


Ignore:
Timestamp:
May 9, 2008, 2:20:08 PM (17 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:
0b05147
Parents:
de293ac
Message:

MinimumRingSize is now an array over all atoms

Each entry in MinimumRingSize defines the maximum bond order possible due to membership in or vicinity to a ring. In order to do create this array, CyclicStructureAnalysis() was expanded with another BFS going over all atoms that are only close to a ring (nearest ring member is search and entri in array set). All other functions such as DepthFirstSearchAnalysis() and FragmentBOSSANOVA(), where the MinimumRingSize is now checked (and not in FragmentMolecule) anymore, are adapted to have an array parameter.

Location:
src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    rde293ac rfc850d  
    720720  int Z;
    721721  int j, axis, count, faktor;
    722   int MinimumRingSize = -1;
     722  int *MinimumRingSize = NULL;
    723723  enum ConfigStatus config_present = absent;
    724724  MoleculeLeafClass *Subgraphs = NULL;
    … …  
    11401140        }
    11411141        delete(Subgraphs);    // we don't need the list here, so free everything
     1142        delete[](MinimumRingSize);
    11421143        Subgraphs = NULL;
    11431144        end = clock();

  • src/molecules.cpp

    rde293ac rfc850d  
    11791179{
    11801180  int No = 0;
    1181   int MinimumRingSize;
     1181  int *MinimumRingSize = NULL;
    11821182  MoleculeLeafClass *Subgraphs = NULL;
    11831183  bond *Binder = first;
    … …  
    11901190    }
    11911191    delete(Subgraphs);
     1192    delete[](MinimumRingSize);
    11921193  }
    11931194  while(Binder->next != last) {
    … …  
    14371438 * \param *out output stream for debugging
    14381439 * \param ReturnStack true - return pointer to atom stack of separable components, false - return NULL
    1439  * \param MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found
     1440 * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found
    14401441 * \return list of each disconnected subgraph as an individual molecule class structure
    14411442 */
    1442 MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int &MinimumRingSize)
     1443MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize)
    14431444{
    14441445  class StackClass<atom *> *AtomStack;
    … …  
    16301631 * \param *out output stream for debugging
    16311632 * \param *BackEdgeStack stack with all back edges found during DFS scan
    1632  * \param MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
     1633 * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
    16331634 * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
    16341635 */
    1635 void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *  BackEdgeStack, int &MinimumRingSize)
     1636void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *  BackEdgeStack, int *&MinimumRingSize)
    16361637{
    16371638  atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CyclicStructureAnalysis: **PredecessorList");
    … …  
    16421643  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
    16431644  bond *Binder = NULL, *BackEdge = NULL;
    1644   int RingSize, NumCycles;
     1645  int RingSize, NumCycles, MinRingSize = -1;
    16451646
    16461647  // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
    … …  
    16501651    ColorList[i] = white;
    16511652  }
     1653  MinimumRingSize = new int[AtomCount];
     1654  for(int i=0;i<AtomCount;i++)
     1655    MinimumRingSize[i] = AtomCount;
     1656
    16521657 
    16531658  *out << Verbose(1) << "Back edge list - ";
    … …  
    16551660 
    16561661  *out << Verbose(1) << "Analysing cycles ... " << endl;
    1657   if ((MinimumRingSize <= 2) || (MinimumRingSize > AtomCount))
    1658     MinimumRingSize = AtomCount;
    16591662  NumCycles = 0;
    16601663  while (!BackEdgeStack->IsEmpty()) {
    … …  
    16701673    //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
    16711674    OtherAtom = NULL;
    1672     while ((Walker != Root) && ((OtherAtom == NULL) || (ShortestPathList[OtherAtom->nr] < MinimumRingSize))) {  // look for Root
     1675    while ((Walker != Root) && ((OtherAtom == NULL) || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]))) {  // look for Root
    16731676      Walker = BFSStack->PopFirst();
    16741677      //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
    … …  
    16841687            ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
    16851688            //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
    1686             if (ShortestPathList[OtherAtom->nr] < MinimumRingSize) { // Check for maximum distance
     1689            if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
    16871690              //*out << Verbose(3) << "Putting OtherAtom into queue." << endl;
    16881691              BFSStack->Push(OtherAtom);
    … …  
    17111714      }
    17121715      *out << Walker->Name << "  with a length of " << RingSize << "." << endl << endl;
    1713       if (RingSize < MinimumRingSize)
    1714         MinimumRingSize = RingSize;
     1716      // walk through all and set MinimumRingSize
     1717      Walker = Root;
     1718      while (Walker != BackEdge->rightatom) {
     1719        Walker = PredecessorList[Walker->nr];
     1720        if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
     1721          MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
     1722      }
     1723      if ((RingSize < MinRingSize) || (MinRingSize == -1))
     1724        MinRingSize = RingSize;
    17151725    } else {
    1716       *out << Verbose(1) << "No ring with length equal to or smaller than " << MinimumRingSize << " found." << endl;
     1726      *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
    17171727    }
    17181728   
    … …  
    17251735    }
    17261736  }
     1737  // go over all atoms
     1738  Root = start;
     1739  while(Root->next != end) {
     1740    Root = Root->next;
     1741   
     1742    if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
     1743      ShortestPathList[Walker->nr] = 0;
     1744      BFSStack->ClearStack();  // start with empty BFS stack
     1745      BFSStack->Push(Walker);
     1746      TouchedStack->Push(Walker);
     1747      //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
     1748      OtherAtom = Walker;
     1749      while ((Walker != Root) && (OtherAtom != NULL)) {  // look for Root
     1750        Walker = BFSStack->PopFirst();
     1751        //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
     1752        for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
     1753          Binder = ListOfBondsPerAtom[Walker->nr][i];
     1754          if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
     1755            OtherAtom = Binder->GetOtherAtom(Walker);
     1756            //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
     1757            if (ColorList[OtherAtom->nr] == white) {
     1758              TouchedStack->Push(OtherAtom);
     1759              ColorList[OtherAtom->nr] = lightgray;
     1760              PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
     1761              ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
     1762              //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
     1763              if (MinimumRingSize[OtherAtom->GetTrueFather()->nr] != AtomCount) { // if the other atom is connected to a ring
     1764                MinimumRingSize[Root->GetTrueFather()->nr] = ShortestPathList[OtherAtom->nr]+MinimumRingSize[OtherAtom->GetTrueFather()->nr];
     1765                OtherAtom = NULL; //break;
     1766                break;
     1767              } else
     1768                BFSStack->Push(OtherAtom);
     1769            } else {
     1770              //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
     1771            }
     1772          } else {
     1773            //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
     1774          }
     1775        }
     1776        ColorList[Walker->nr] = black;
     1777        //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
     1778      }
     1779
     1780      // now clean the lists
     1781      while (!TouchedStack->IsEmpty()){
     1782        Walker = TouchedStack->PopFirst();
     1783        PredecessorList[Walker->nr] = NULL;
     1784        ShortestPathList[Walker->nr] = -1;
     1785        ColorList[Walker->nr] = white;
     1786      }
     1787    }
     1788    *out << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl;
     1789  }
     1790 
    17271791       
    1728 
    1729   if (MinimumRingSize != -1)
    1730     *out << Verbose(1) << "Minimum ring size is " << MinimumRingSize << ", over " << NumCycles << " cycles total." << endl;
     1792  if (MinRingSize != -1)
     1793    *out << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl;
    17311794  else
    17321795    *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
    … …  
    21482211  bool status = false;
    21492212  ifstream InputFile;
    2150   stringstream line;
     2213
     2214  // initialize mask list
     2215  for(int i=0;i<AtomCount;i++)
     2216    AtomMask[i] = false;
    21512217 
    21522218  if (Order < 0) { // adaptive increase of BondOrder per site
    2153     for(int i=0;i<AtomCount;i++)
    2154       AtomMask[i] = false;
     2219    if (AtomMask[AtomCount] == true)  // break after one step
     2220      return false;
    21552221    // parse the EnergyPerFragment file
    2156     line.str(path);
    2157     line << "/" << FRAGMENTPREFIX << ENERGYPERFRAGMENT;
    2158     InputFile.open(line.str().c_str(), ios::in);
    2159     if (InputFile != NULL) {
    2160       char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
     2222    char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
     2223    sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
     2224    InputFile.open(buffer, ios::in);
     2225    if ((InputFile != NULL) && (GlobalKeySetList != NULL)) {
     2226      // transmorph graph keyset list into indexed KeySetList
     2227      map<int,KeySet> IndexKeySetList;
     2228      for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
     2229        IndexKeySetList.insert( pair<int,KeySet>(runner->second.first,runner->first) );
     2230      }
    21612231      int lines = 0;
    21622232      // count the number of lines, i.e. the number of fragments
     2233      InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
     2234      InputFile.getline(buffer, MAXSTRINGSIZE);
    21632235      while(!InputFile.eof()) {
    21642236        InputFile.getline(buffer, MAXSTRINGSIZE);
    21652237        lines++;
    21662238      }
     2239      *out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
    21672240      InputFile.clear();
    21682241      InputFile.seekg(ios::beg);
    2169       map<double, int> FragmentEnergies;
    2170       int No;
    2171       int Value;
     2242      map<int, pair<double,int> > AdaptiveCriteriaList;  // (Root No., (Value, Order)) !
     2243      int No, FragOrder;
     2244      double Value;
    21722245      // each line represents a fragment root (Atom::nr) id and its energy contribution
     2246      InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
     2247      InputFile.getline(buffer, MAXSTRINGSIZE);
    21732248      while(!InputFile.eof()) {
    21742249        InputFile.getline(buffer, MAXSTRINGSIZE);
    2175         InputFile >> No;
    2176         InputFile >> Value;
    2177         FragmentEnergies.insert( make_pair(Value, No) );
    2178       }
    2179       for(map<double, int>::iterator runner = FragmentEnergies.lower_bound(pow(1.,Order)); runner != FragmentEnergies.begin(); runner++) {
    2180         // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
    2181         AtomMask[(*runner).second] = true;
     2250        if (strlen(buffer) > 2) {
     2251          //*out << Verbose(2) << "Scanning: " << buffer;
     2252          stringstream line(buffer);
     2253          line >> FragOrder;
     2254          line >> ws >> No;
     2255          line >> ws >> Value; // skip time entry
     2256          line >> ws >> Value;
     2257          No -= 1;  // indices start at 1 in file, not 0
     2258          //*out << Verbose(2) << " - yields (" << No << "," << Value << ")" << endl;
     2259
     2260          // clean the list of those entries that have been superceded by higher order terms already
     2261          map<int,KeySet>::iterator marker = IndexKeySetList.find(No);    // find keyset to Frag No.
     2262          if (marker != IndexKeySetList.end()) {  // if found
     2263            // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
     2264            pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( Value, Order) ));
     2265            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
     2266            if (!InsertedElement.second) { // this root is already present
     2267              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
     2268                //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
     2269                {  // if value is smaller, update value and order
     2270                (*PresentItem).second.first = Value;
     2271                (*PresentItem).second.second = FragOrder;
     2272                *out << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
     2273              }
     2274            } else {
     2275              *out << Verbose(2) << "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
     2276            }
     2277          } else {
     2278            *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
     2279          }
     2280        }
     2281      }
     2282      // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
     2283      map<double, pair<int,int> > FinalRootCandidates;
     2284      *out << Verbose(1) << "Root candidate list is: " << endl;
     2285      for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList.begin(); runner != AdaptiveCriteriaList.end(); runner++) {
     2286        Walker = FindAtom((*runner).first);
     2287        if (Walker != NULL) {
     2288          if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
     2289            *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
     2290            FinalRootCandidates.insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
     2291          }
     2292        } else {
     2293          cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
     2294        }
     2295      }
     2296      // pick the ones still below threshold and mark as to be adaptively updated
     2297      for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
     2298        No = (*runner).second.first;
     2299        *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
     2300        AtomMask[No] = true;
     2301        status = true;
    21822302      }
    21832303      // close and done
    21842304      InputFile.close();
    21852305      InputFile.clear();
    2186       Free((void **)&buffer, "molecule::CheckOrderAtSite: *buffer");
    21872306    } else {
    2188       cerr << "Unable to parse " << FRAGMENTPREFIX << ENERGYPERFRAGMENT << " file." << endl;
    2189       return false;
    2190     }
     2307      cerr << "Unable to parse " << buffer << " file, incrementing all." << endl;
     2308      while (Walker->next != end) {
     2309        Walker = Walker->next;
     2310    #ifdef ADDHYDROGEN
     2311        if (Walker->type->Z != 1) // skip hydrogen
     2312    #endif
     2313        {
     2314          AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
     2315          status = true;
     2316        }
     2317      }
     2318    }
     2319    Free((void **)&buffer, "molecule::CheckOrderAtSite: *buffer");
    21912320    // pick a given number of highest values and set AtomMask
    2192   } else if (Order > 0) { // global increase of Bond Order
     2321  } else { // global increase of Bond Order
    21932322    while (Walker->next != end) {
    21942323      Walker = Walker->next;
    … …  
    21982327      {
    21992328        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
    2200         if (Walker->AdaptiveOrder < Order)
     2329        if ((Order != 0) && (Walker->AdaptiveOrder < Order))
    22012330          status = true;
    22022331      }
    22032332    }
     2333    if ((Order == 0) && (AtomMask[AtomCount] == true))  // single stepping, just check
     2334      status = false;
     2335     
    22042336    if (!status)
    22052337      *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
    2206   } else {  // single stepping, just check
    2207     if (AtomMask[AtomCount] != true)  // if true, we would've done a step already
    2208       for(int i=0;i<=AtomCount;i++)
    2209         AtomMask[i] = true;
    2210     else
    2211       status = false;
    2212   }
     2338  }
     2339 
     2340  // print atom mask for debugging
     2341  *out << "              ";
     2342  for(int i=0;i<AtomCount;i++)
     2343    *out << (i % 10);
     2344  *out << endl << "Atom mask is: ";
     2345  for(int i=0;i<AtomCount;i++)
     2346    *out << (AtomMask[i] ? "t" : "f");
     2347  *out << endl;
    22132348 
    22142349  return status;
    … …  
    22682403        MoleculeListClass *BondFragments = NULL;
    22692404  int *SortIndex = NULL;
    2270   int MinimumRingSize;
     2405  int *MinimumRingSize = NULL;
    22712406  int FragmentCounter;
    22722407  MoleculeLeafClass *MolecularWalker = NULL;
    … …  
    23152450 
    23162451  // =================================== Begin of FRAGMENTATION ===============================
    2317   // check cyclic lengths
    2318   if ((MinimumRingSize != -1) && (Order >= MinimumRingSize)) {
    2319     *out << Verbose(0) << "Bond order " << Order << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
    2320   } else {
    2321     // ===== 6a. assign each keyset to its respective subgraph =====
    2322     Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), false);
    2323 
    2324     KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
    2325     AtomMask = new bool[AtomCount+1];
    2326     for (int i=0;i<AtomCount;i++) // need to set to false to recognize in single-stepping (one global step set them all to true)
    2327       AtomMask[i] = false;
    2328     while (CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order)) {
    2329       AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    2330       // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    2331       Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, Order, (FragmentCounter = 0));
    2332 
    2333       // ===== 7. fill the bond fragment list =====
    2334       FragmentCounter = 0;
    2335       MolecularWalker = Subgraphs;
    2336       while (MolecularWalker->next != NULL) {
    2337         MolecularWalker = MolecularWalker->next;
    2338         *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
    2339         // output ListOfBondsPerAtom for debugging
    2340         MolecularWalker->Leaf->OutputListOfBonds(out);
    2341         if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
    2342        
    2343           // call BOSSANOVA method
    2344           *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
    2345           MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter]);
    2346         } else {
    2347           cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
    2348         }
    2349         FragmentCounter++;  // next fragment list
    2350       }
    2351     }
    2352     delete[](RootStack);
    2353     delete[](AtomMask);
    2354   }
     2452  // ===== 6a. assign each keyset to its respective subgraph =====
     2453  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), false);
     2454
     2455  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
     2456  AtomMask = new bool[AtomCount+1];
     2457  while (CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, configuration->configpath)) {
     2458    AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
     2459    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
     2460    Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
     2461
     2462    // ===== 7. fill the bond fragment list =====
     2463    FragmentCounter = 0;
     2464    MolecularWalker = Subgraphs;
     2465    while (MolecularWalker->next != NULL) {
     2466      MolecularWalker = MolecularWalker->next;
     2467      *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
     2468      // output ListOfBondsPerAtom for debugging
     2469      MolecularWalker->Leaf->OutputListOfBonds(out);
     2470      if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
     2471     
     2472        // call BOSSANOVA method
     2473        *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
     2474        MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
     2475      } else {
     2476        cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
     2477      }
     2478      FragmentCounter++;  // next fragment list
     2479    }
     2480  }
     2481  delete[](RootStack);
     2482  delete[](AtomMask);
    23552483  delete(ParsedFragmentList);
     2484  delete[](MinimumRingSize);
    23562485 
    23572486  // free the index lookup list
    … …  
    36273756 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
    36283757 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
     3758 * \param *MinimumRingSize minimum ring size for each atom (molecule::Atomcount)
    36293759 * \return pointer to Graph list
    36303760 */
    3631 void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack)
     3761void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
    36323762{
    36333763  Graph ***FragmentLowerOrdersList = NULL;
    … …  
    36753805    RootKeyNr = RootStack.front();
    36763806    RootStack.pop_front();
    3677     // increase adaptive order by one
    36783807    Walker = FindAtom(RootKeyNr);
    3679     Walker->GetTrueFather()->AdaptiveOrder++;
    3680     Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
    3681 
    3682     // initialise Order-dependent entries of UniqueFragments structure
    3683     FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
    3684     FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
    3685     for (int i=0;i<Order;i++) {
    3686       FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
    3687       FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
    3688       FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
    3689       FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
    3690       FragmentSearch.BondsPerSPCount[i] = 0;
    3691     } 
    3692 
    3693     // allocate memory for all lower level orders in this 1D-array of ptrs
    3694     NumLevels = 1 << (Order-1); // (int)pow(2,Order);
    3695     FragmentLowerOrdersList[RootNr] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
    3696  
    3697     // create top order where nothing is reduced
    3698     *out << Verbose(0) << "==============================================================================================================" << endl;
    3699     *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr-1) << " Roots remaining." << endl;
    3700 
    3701     // Create list of Graphs of current Bond Order (i.e. F_{ij})
    3702     FragmentLowerOrdersList[RootNr][0] =  new Graph;
    3703     FragmentSearch.TEFactor = 1.;
    3704     FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0];      // set to insertion graph
    3705     FragmentSearch.Root = Walker;
    3706     NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
    3707     *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    3708     NumMolecules = 0;
     3808    // check cyclic lengths
     3809    if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 >= MinimumRingSize[Walker->GetTrueFather()->nr])) {
     3810      *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;
     3811    } else {
     3812      // increase adaptive order by one
     3813      Walker->GetTrueFather()->AdaptiveOrder++;
     3814      Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
     3815 
     3816      // initialise Order-dependent entries of UniqueFragments structure
     3817      FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
     3818      FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
     3819      for (int i=0;i<Order;i++) {
     3820        FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
     3821        FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
     3822        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
     3823        FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
     3824        FragmentSearch.BondsPerSPCount[i] = 0;
     3825      } 
     3826 
     3827      // allocate memory for all lower level orders in this 1D-array of ptrs
     3828      NumLevels = 1 << (Order-1); // (int)pow(2,Order);
     3829      FragmentLowerOrdersList[RootNr] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
    37093830   
    3710     if ((NumLevels >> 1) > 0) {
    3711       // create lower order fragments
    3712       *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
    3713       Order = Walker->AdaptiveOrder;
    3714       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)
    3715         // step down to next order at (virtual) boundary of powers of 2 in array
    3716         while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
    3717           Order--;
    3718         *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
    3719         for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
    3720           int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
    3721           *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
    3722           *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
    3723    
    3724           // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
    3725           //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
    3726           //NumMolecules = 0;
    3727           FragmentLowerOrdersList[RootNr][dest] = new Graph;
    3728           for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
    3729             for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    3730               Graph TempFragmentList;
    3731               FragmentSearch.TEFactor = -(*runner).second.second;
    3732               FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
    3733               FragmentSearch.Root = FindAtom(*sprinter);
    3734               NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
    3735               // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
    3736               *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
    3737               InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
     3831      // create top order where nothing is reduced
     3832      *out << Verbose(0) << "==============================================================================================================" << endl;
     3833      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr-1) << " Roots remaining." << endl;
     3834 
     3835      // Create list of Graphs of current Bond Order (i.e. F_{ij})
     3836      FragmentLowerOrdersList[RootNr][0] =  new Graph;
     3837      FragmentSearch.TEFactor = 1.;
     3838      FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0];      // set to insertion graph
     3839      FragmentSearch.Root = Walker;
     3840      NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
     3841      *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     3842      NumMolecules = 0;
     3843     
     3844      if ((NumLevels >> 1) > 0) {
     3845        // create lower order fragments
     3846        *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
     3847        Order = Walker->AdaptiveOrder;
     3848        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)
     3849          // step down to next order at (virtual) boundary of powers of 2 in array
     3850          while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
     3851            Order--;
     3852          *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
     3853          for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
     3854            int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
     3855            *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
     3856            *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
     3857     
     3858            // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
     3859            //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
     3860            //NumMolecules = 0;
     3861            FragmentLowerOrdersList[RootNr][dest] = new Graph;
     3862            for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
     3863              for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
     3864                Graph TempFragmentList;
     3865                FragmentSearch.TEFactor = -(*runner).second.second;
     3866                FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
     3867                FragmentSearch.Root = FindAtom(*sprinter);
     3868                NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
     3869                // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
     3870                *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
     3871                InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
     3872              }
    37383873            }
     3874            *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
    37393875          }
    3740           *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
    37413876        }
    37423877      }
    3743     }
    3744     // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
    3745     //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
    3746     TotalNumMolecules += NumMoleculesOfOrder[RootNr];
    3747     *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    3748     RootStack.push_back(RootKeyNr); // put back on stack
    3749     RootNr++;
    3750 
    3751     // free Order-dependent entries of UniqueFragments structure for next loop cycle
    3752     Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
    3753     for (int i=0;i<Order;i++) {
    3754       delete(FragmentSearch.BondsPerSPList[2*i]);
    3755       delete(FragmentSearch.BondsPerSPList[2*i+1]);
    3756     } 
    3757     Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
     3878      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
     3879      //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
     3880      TotalNumMolecules += NumMoleculesOfOrder[RootNr];
     3881      *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     3882      RootStack.push_back(RootKeyNr); // put back on stack
     3883      RootNr++;
     3884 
     3885      // free Order-dependent entries of UniqueFragments structure for next loop cycle
     3886      Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
     3887      for (int i=0;i<Order;i++) {
     3888        delete(FragmentSearch.BondsPerSPList[2*i]);
     3889        delete(FragmentSearch.BondsPerSPList[2*i+1]);
     3890      } 
     3891      Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
     3892    }
    37583893  }
    37593894  *out << Verbose(0) << "==============================================================================================================" << endl;

  • src/molecules.hpp

    rde293ac rfc850d  
    326326 
    327327  // Graph analysis
    328   MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int &MinimumRingSize);
    329   void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int &MinimumRingSize);
     328  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize);
     329  void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
    330330  bond * FindNextUnused(atom *vertex);
    331331  void SetNextComponentNumber(atom *vertex, int nr);
    … …  
    353353  void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
    354354  /// -# BOSSANOVA
    355   void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack);
     355  void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
    356356        int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
    357357  bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
    … …  
    413413  bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
    414414  bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
    415   bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int Order, int &FragmentCounter);
     415  bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
    416416  bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
    417417  bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, int &FragmentCounter, int GlobalAtomCount, bool &FreeList);
Note: See TracChangeset for help on using the changeset viewer.