Changeset a3ff7b2 for molecuilder


Ignore:
Timestamp:
May 8, 2008, 5:31:04 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
310b25
Parents:
461e93
Message:

molecule::CheckOrderAtSite() now interprets negative Orders as adaptive increase by parsing EnergyPerFragment.da
t

positive Order is global increase till all sites have at least this order, Order 0 means single global increase step, negative Order is the exponent in 10{-Order} to give the threshold value: ENERGYPERFRAGMENT is scanned in
to a map (sorted by values) and all fragments still above the threshold are taken into AtomMask for increase. Jo
iner has to write this file with (Root id of Fragment, energy contribution) pairs.

Location:
molecuilder/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/defs.hpp

    r461e93 ra3ff7b2  
    4646#define TEFACTORSFILE "TE-Factors.dat"
    4747#define FORCESFILE "Forces-Factors.dat"
    48 #define ORDERATSITEFILE "OrderAtSite.dat"
     48#define ORDERATSITEFILE "OrderAtSite.dat"
     49#define ENERGYPERFRAGMENT "EnergyPerFragment.dat"
    4950#define FRAGMENTPREFIX "BondFragment"
    5051#define STANDARDCONFIG "unknown.conf"

  • molecuilder/src/moleculelist.cpp

    r461e93 ra3ff7b2  
    408408 * \param *out output stream for debugging
    409409 * \param *&RootStack stack to be filled
     410 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site
    410411 * \param Order desired bond order for all sites
    411412 * \param &FragmentCounter counts through the fragments in this MoleculeLeafClass
    412413 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
    413414 */
    414 bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, int Order, int &FragmentCounter)
    415 {
    416   atom *Walker = NULL;
     415bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int Order, int &FragmentCounter)
     416{
     417  atom *Walker = NULL, *Father = NULL;
    417418
    418419  if (RootStack != NULL) {
    … …  
    426427        while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
    427428          Walker = Walker->next;
     429          Father = Walker->GetTrueFather();
     430          if (AtomMask[Father->nr]) // apply mask
    428431      #ifdef ADDHYDROGEN
    429           if (Walker->type->Z != 1) // skip hydrogen
     432            if (Walker->type->Z != 1) // skip hydrogen
    430433      #endif
    431             if (Walker->GetTrueFather()->AdaptiveOrder < Order) // only if Order is still greater
    432               RootStack[FragmentCounter].push_front(Walker->nr);
     434              if (Father->AdaptiveOrder < Order) // only if Order is still greater
     435                RootStack[FragmentCounter].push_front(Walker->nr);
    433436        }
    434437        if (next != NULL)
    435           next->FillRootStackForSubgraphs(out, RootStack, Order, ++FragmentCounter);
     438          next->FillRootStackForSubgraphs(out, RootStack, AtomMask, Order, ++FragmentCounter);
    436439      }  else {
    437440        *out << Verbose(1) << "Rootstack[" << FragmentCounter  << "] is NULL." << endl;

  • molecuilder/src/molecules.cpp

    r461e93 ra3ff7b2  
    20172017    *out << Verbose(1) << "done." << endl;
    20182018  } else {
    2019     *out << Verbose(1) << "failed to open file" << line << "." << endl;
     2019    *out << Verbose(1) << "failed to open " << line << "." << endl;
    20202020    status = false;
    20212021  }
    … …  
    21372137/** Checks whether the OrderAtSite is still bewloe \a Order at some site.
    21382138 * \param *out output stream for debugging
    2139  * \param &FragmentationToDo return boolean
    2140  * \param Order desired Order
     2139 * \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
     2140 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
     2141 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
     2142 * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
    21412143 * \return true - needs further fragmentation, false - does not need fragmentation
    21422144 */
    2143 bool molecule::CheckOrderAtSite(ofstream *out, int Order)
    2144 {
    2145   atom *Walker = NULL;
     2145bool molecule::CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, char *path)
     2146{
     2147  atom *Walker = start;
    21462148  bool status = false;
    2147  
    2148   Walker = start;
    2149   while (Walker->next != end) {
    2150     Walker = Walker->next;
    2151 #ifdef ADDHYDROGEN
    2152     if (Walker->type->Z != 1) // skip hydrogen
    2153 #endif
    2154       if (Walker->AdaptiveOrder < Order)
    2155         status = true;
    2156   }
    2157   if (!status)
    2158     *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
     2149  ifstream InputFile;
     2150  stringstream line;
     2151 
     2152  if (Order < 0) { // adaptive increase of BondOrder per site
     2153    for(int i=0;i<AtomCount;i++)
     2154      AtomMask[i] = false;
     2155    // 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");
     2161      int lines = 0;
     2162      // count the number of lines, i.e. the number of fragments
     2163      while(!InputFile.eof()) {
     2164        InputFile.getline(buffer, MAXSTRINGSIZE);
     2165        lines++;
     2166      }
     2167      InputFile.clear();
     2168      InputFile.seekg(ios::beg);
     2169      map<double, int> FragmentEnergies;
     2170      int No;
     2171      int Value;
     2172      // each line represents a fragment root (Atom::nr) id and its energy contribution
     2173      while(!InputFile.eof()) {
     2174        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;
     2182      }
     2183      // close and done
     2184      InputFile.close();
     2185      InputFile.clear();
     2186      Free((void **)&buffer, "molecule::CheckOrderAtSite: *buffer");
     2187    } else {
     2188      cerr << "Unable to parse " << FRAGMENTPREFIX << ENERGYPERFRAGMENT << " file." << endl;
     2189      return false;
     2190    }
     2191    // pick a given number of highest values and set AtomMask
     2192  } else if (Order > 0) { // global increase of Bond Order
     2193    while (Walker->next != end) {
     2194      Walker = Walker->next;
     2195  #ifdef ADDHYDROGEN
     2196      if (Walker->type->Z != 1) // skip hydrogen
     2197  #endif
     2198      {
     2199        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
     2200        if (Walker->AdaptiveOrder < Order)
     2201          status = true;
     2202      }
     2203    }
     2204    if (!status)
     2205      *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  }
    21592213 
    21602214  return status;
    … …  
    22262280  atom **ListOfAtoms = NULL;
    22272281  atom ***ListOfLocalAtoms = NULL;
     2282  bool *AtomMask = NULL;
    22282283 
    22292284  *out << endl;
    … …  
    22662321    // ===== 6a. assign each keyset to its respective subgraph =====
    22672322    Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), false);
    2268     delete(ParsedFragmentList);
    22692323
    22702324    KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
    2271     while (CheckOrderAtSite(out, Order)) {
     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()
    22722330      // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    2273       Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, Order, (FragmentCounter = 0));
     2331      Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, Order, (FragmentCounter = 0));
    22742332
    22752333      // ===== 7. fill the bond fragment list =====
    … …  
    22932351    }
    22942352    delete[](RootStack);
    2295   }
     2353    delete[](AtomMask);
     2354  }
     2355  delete(ParsedFragmentList);
    22962356 
    22972357  // free the index lookup list

  • molecuilder/src/molecules.hpp

    r461e93 ra3ff7b2  
    3737#define KeyStack deque<int>
    3838#define KeySet set<int>
    39 #define Graph map<KeySet, pair<int, double>, KeyCompare >
    40 #define GraphPair pair<KeySet, pair<int, double> >
     39#define NumberValuePair pair<int, double>
     40#define Graph map<KeySet, NumberValuePair, KeyCompare >
     41#define GraphPair pair<KeySet, NumberValuePair >
    4142#define KeySetTestPair pair<KeySet::iterator, bool>
    4243#define GraphTestPair pair<Graph::iterator, bool>
    … …  
    340341  /// Fragment molecule by two different approaches:
    341342  void FragmentMolecule(ofstream *out, int Order, config *configuration);
    342   bool CheckOrderAtSite(ofstream *out, int Order);
     343  bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, char *path = NULL);
    343344  bool StoreAdjacencyToFile(ofstream *out, char *path);
    344345  bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
    … …  
    412413  bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
    413414  bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
    414   bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, int Order, int &FragmentCounter);
     415  bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int Order, int &FragmentCounter);
    415416  bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
    416417  bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, int &FragmentCounter, int GlobalAtomCount, bool &FreeList);
Note: See TracChangeset for help on using the changeset viewer.