| [14de469] | 1 | /** \file MoleculeListClass.cpp
 | 
|---|
 | 2 |  * 
 | 
|---|
 | 3 |  * Function implementations for the class MoleculeListClass.
 | 
|---|
 | 4 |  * 
 | 
|---|
 | 5 |  */
 | 
|---|
 | 6 | 
 | 
|---|
 | 7 | #include "molecules.hpp"
 | 
|---|
 | 8 | 
 | 
|---|
 | 9 | /*********************************** Functions for class MoleculeListClass *************************/
 | 
|---|
 | 10 | 
 | 
|---|
 | 11 | /** Constructor for MoleculeListClass.
 | 
|---|
 | 12 |  */
 | 
|---|
 | 13 | MoleculeListClass::MoleculeListClass()
 | 
|---|
 | 14 | {
 | 
|---|
 | 15 | };
 | 
|---|
 | 16 | 
 | 
|---|
 | 17 | /** constructor for MoleculeListClass.
 | 
|---|
 | 18 |  * \param NumMolecules number of molecules to allocate for
 | 
|---|
 | 19 |  * \param NumAtoms number of atoms to allocate for
 | 
|---|
 | 20 |  */
 | 
|---|
 | 21 | MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
 | 
|---|
 | 22 | {
 | 
|---|
 | 23 |   ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
 | 
|---|
| [7f3b9d] | 24 |   for (int i=NumMolecules;i--;)
 | 
|---|
| [14de469] | 25 |     ListOfMolecules[i] = NULL;
 | 
|---|
 | 26 |   NumberOfMolecules = NumMolecules;
 | 
|---|
 | 27 |   NumberOfTopAtoms = NumAtoms;
 | 
|---|
 | 28 | };
 | 
|---|
 | 29 | 
 | 
|---|
 | 30 | 
 | 
|---|
 | 31 | /** Destructor for MoleculeListClass.
 | 
|---|
 | 32 |  */
 | 
|---|
 | 33 | MoleculeListClass::~MoleculeListClass()
 | 
|---|
 | 34 | {
 | 
|---|
| [db942e] | 35 |   cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
 | 
|---|
| [7f3b9d] | 36 |   for (int i=NumberOfMolecules;i--;) {
 | 
|---|
| [14de469] | 37 |     if (ListOfMolecules[i] != NULL) { // if NULL don't free
 | 
|---|
 | 38 |       cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
 | 
|---|
 | 39 |       delete(ListOfMolecules[i]);
 | 
|---|
 | 40 |       ListOfMolecules[i] = NULL;
 | 
|---|
 | 41 |     }
 | 
|---|
 | 42 |   }
 | 
|---|
 | 43 |   cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
 | 
|---|
 | 44 |   Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
 | 
|---|
 | 45 | };
 | 
|---|
 | 46 | 
 | 
|---|
| [db942e] | 47 | /** Compare whether two molecules are equal.
 | 
|---|
 | 48 |  * \param *a molecule one
 | 
|---|
 | 49 |  * \param *n molecule two
 | 
|---|
 | 50 |  * \return lexical value (-1, 0, +1)
 | 
|---|
 | 51 |  */
 | 
|---|
| [14de469] | 52 | int MolCompare(const void *a, const void *b)
 | 
|---|
 | 53 | {
 | 
|---|
 | 54 |   int *aList = NULL, *bList = NULL;
 | 
|---|
 | 55 |   int Count, Counter, aCounter, bCounter;
 | 
|---|
 | 56 |   int flag;
 | 
|---|
 | 57 |   atom *aWalker = NULL;
 | 
|---|
 | 58 |   atom *bWalker = NULL;
 | 
|---|
 | 59 |   
 | 
|---|
 | 60 |   // sort each atom list and put the numbers into a list, then go through
 | 
|---|
 | 61 |   //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
 | 
|---|
 | 62 |   if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
 | 
|---|
 | 63 |     return -1;
 | 
|---|
 | 64 |   } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
 | 
|---|
 | 65 |     return +1;
 | 
|---|
 | 66 |     else {
 | 
|---|
 | 67 |       Count = (**(molecule **)a).AtomCount;
 | 
|---|
| [7f3b9d] | 68 |       aList = new int[Count];
 | 
|---|
 | 69 |       bList = new int[Count];
 | 
|---|
| [14de469] | 70 |   
 | 
|---|
 | 71 |       // fill the lists
 | 
|---|
 | 72 |       aWalker = (**(molecule **)a).start;
 | 
|---|
 | 73 |       bWalker = (**(molecule **)b).start;
 | 
|---|
 | 74 |       Counter = 0;
 | 
|---|
 | 75 |       aCounter = 0;
 | 
|---|
 | 76 |       bCounter = 0;
 | 
|---|
 | 77 |       while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
 | 
|---|
 | 78 |         aWalker = aWalker->next;
 | 
|---|
 | 79 |         bWalker = bWalker->next;
 | 
|---|
 | 80 |         if (aWalker->GetTrueFather() == NULL)
 | 
|---|
 | 81 |           aList[Counter] = Count + (aCounter++);
 | 
|---|
 | 82 |         else
 | 
|---|
 | 83 |           aList[Counter] = aWalker->GetTrueFather()->nr;
 | 
|---|
 | 84 |         if (bWalker->GetTrueFather() == NULL)
 | 
|---|
 | 85 |           bList[Counter] = Count + (bCounter++);
 | 
|---|
 | 86 |         else
 | 
|---|
 | 87 |           bList[Counter] = bWalker->GetTrueFather()->nr;
 | 
|---|
 | 88 |         Counter++;
 | 
|---|
 | 89 |       }
 | 
|---|
 | 90 |       // check if AtomCount was for real
 | 
|---|
 | 91 |       flag = 0;
 | 
|---|
 | 92 |       if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
 | 
|---|
 | 93 |         flag = -1;
 | 
|---|
 | 94 |       } else {
 | 
|---|
 | 95 |         if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
 | 
|---|
 | 96 |           flag = 1;
 | 
|---|
 | 97 |       }
 | 
|---|
 | 98 |       if (flag == 0) {
 | 
|---|
 | 99 |         // sort the lists
 | 
|---|
 | 100 |         gsl_heapsort(aList,Count, sizeof(int), CompareDoubles); 
 | 
|---|
 | 101 |         gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
 | 
|---|
 | 102 |         // compare the lists 
 | 
|---|
 | 103 |         
 | 
|---|
 | 104 |         flag = 0;
 | 
|---|
 | 105 |         for(int i=0;i<Count;i++) {
 | 
|---|
 | 106 |           if (aList[i] < bList[i]) {
 | 
|---|
 | 107 |             flag = -1;
 | 
|---|
 | 108 |           } else {
 | 
|---|
 | 109 |             if (aList[i] > bList[i])
 | 
|---|
 | 110 |               flag = 1;
 | 
|---|
 | 111 |           }
 | 
|---|
 | 112 |           if (flag != 0)
 | 
|---|
 | 113 |             break;
 | 
|---|
 | 114 |         }
 | 
|---|
 | 115 |       }
 | 
|---|
| [7f3b9d] | 116 |       delete[](aList);
 | 
|---|
 | 117 |       delete[](bList);
 | 
|---|
| [14de469] | 118 |       return flag;
 | 
|---|
 | 119 |     }
 | 
|---|
 | 120 |   }
 | 
|---|
 | 121 |   return  -1;
 | 
|---|
 | 122 | };
 | 
|---|
 | 123 | 
 | 
|---|
 | 124 | /** Simple output of the pointers in ListOfMolecules.
 | 
|---|
 | 125 |  * \param *out output stream
 | 
|---|
 | 126 |  */
 | 
|---|
 | 127 | void MoleculeListClass::Output(ofstream *out)
 | 
|---|
 | 128 | {
 | 
|---|
 | 129 |   *out<< Verbose(1) << "MoleculeList: ";
 | 
|---|
 | 130 |   for (int i=0;i<NumberOfMolecules;i++)
 | 
|---|
 | 131 |     *out << ListOfMolecules[i] << "\t";
 | 
|---|
 | 132 |   *out << endl;
 | 
|---|
 | 133 | };
 | 
|---|
 | 134 | 
 | 
|---|
| [2459b1] | 135 | 
 | 
|---|
 | 136 | /** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
 | 
|---|
 | 137 |  * \param *out output stream for debugging
 | 
|---|
 | 138 |  * \param *path path to file
 | 
|---|
 | 139 |  * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
 | 
|---|
 | 140 |  * \return true - file written successfully, false - writing failed
 | 
|---|
 | 141 |  */
 | 
|---|
 | 142 | bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex)
 | 
|---|
 | 143 | {
 | 
|---|
 | 144 |   bool status = true;
 | 
|---|
 | 145 |   ofstream ForcesFile;
 | 
|---|
 | 146 |   stringstream line;
 | 
|---|
 | 147 |   atom *Walker = NULL;
 | 
|---|
 | 148 |   element *runner = NULL;
 | 
|---|
 | 149 | 
 | 
|---|
 | 150 |   // open file for the force factors
 | 
|---|
 | 151 |   *out << Verbose(1) << "Saving  force factors ... ";
 | 
|---|
 | 152 |   line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
 | 
|---|
 | 153 |   ForcesFile.open(line.str().c_str(), ios::out);
 | 
|---|
 | 154 |   if (ForcesFile != NULL) {
 | 
|---|
 | 155 |     //cout << Verbose(1) << "Final AtomicForcesList: ";
 | 
|---|
 | 156 |     //output << prefix << "Forces" << endl;
 | 
|---|
 | 157 |     for(int j=0;j<NumberOfMolecules;j++) {
 | 
|---|
 | 158 |       //if (TEList[j] != 0) {
 | 
|---|
 | 159 |       runner = ListOfMolecules[j]->elemente->start;
 | 
|---|
 | 160 |       while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
 | 
|---|
 | 161 |         runner = runner->next;
 | 
|---|
 | 162 |         if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
 | 
|---|
 | 163 |           Walker = ListOfMolecules[j]->start;
 | 
|---|
 | 164 |           while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
 | 
|---|
 | 165 |             Walker = Walker->next;
 | 
|---|
 | 166 |             if (Walker->type->Z == runner->Z) {
 | 
|---|
 | 167 |               if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
 | 
|---|
 | 168 |                 //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
 | 
|---|
 | 169 |                 ForcesFile <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
 | 
|---|
 | 170 |                 } else  // otherwise a -1 to indicate an added saturation hydrogen
 | 
|---|
 | 171 |                   ForcesFile << "-1\t";
 | 
|---|
 | 172 |               }
 | 
|---|
 | 173 |             }
 | 
|---|
 | 174 |           }
 | 
|---|
 | 175 |       }
 | 
|---|
 | 176 |       ForcesFile << endl;
 | 
|---|
 | 177 |     }
 | 
|---|
 | 178 |     ForcesFile.close();
 | 
|---|
 | 179 |     *out << Verbose(1) << "done." << endl;
 | 
|---|
 | 180 |   } else {
 | 
|---|
 | 181 |     status = false;
 | 
|---|
 | 182 |     *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
 | 
|---|
 | 183 |   }
 | 
|---|
 | 184 |   ForcesFile.close();
 | 
|---|
 | 185 |   
 | 
|---|
 | 186 |   return status;
 | 
|---|
 | 187 | };
 | 
|---|
 | 188 | 
 | 
|---|
| [14de469] | 189 | /** Writes a config file for each molecule in the given \a **FragmentList.
 | 
|---|
| [db942e] | 190 |  * \param *out output stream for debugging
 | 
|---|
| [14de469] | 191 |  * \param *configuration standard configuration to attach atoms in fragment molecule to.
 | 
|---|
 | 192 |  * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
 | 
|---|
 | 193 |  * \return true - success (each file was written), false - something went wrong.
 | 
|---|
 | 194 |  */
 | 
|---|
| [db942e] | 195 | bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
 | 
|---|
| [14de469] | 196 | {
 | 
|---|
 | 197 |   ofstream outputFragment;
 | 
|---|
| [c750cc] | 198 |   char FragmentName[MAXSTRINGSIZE];
 | 
|---|
 | 199 |   char PathBackup[MAXSTRINGSIZE];
 | 
|---|
| [14de469] | 200 |   bool result = true;
 | 
|---|
 | 201 |   bool intermediateResult = true;
 | 
|---|
 | 202 |   atom *Walker = NULL;
 | 
|---|
 | 203 |   vector BoxDimension;
 | 
|---|
| [1e9384] | 204 |   char *FragmentNumber = NULL;
 | 
|---|
| [5e4058] | 205 |   char *path = NULL;
 | 
|---|
| [14de469] | 206 |   int FragmentCounter = 0;
 | 
|---|
| [a8aac8d] | 207 |   ofstream output;
 | 
|---|
| [14de469] | 208 |   
 | 
|---|
 | 209 |   // store the fragments as config and as xyz
 | 
|---|
 | 210 |   for(int i=0;i<NumberOfMolecules;i++) {
 | 
|---|
| [db942e] | 211 |     // save default path as it is changed for each fragment
 | 
|---|
| [5e4058] | 212 |     path = configuration->GetDefaultPath();
 | 
|---|
 | 213 |     if (path != NULL)
 | 
|---|
 | 214 |       strcpy(PathBackup, path);
 | 
|---|
 | 215 |     else
 | 
|---|
 | 216 |       cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
 | 
|---|
| [db942e] | 217 | 
 | 
|---|
 | 218 |     // correct periodic
 | 
|---|
 | 219 |     ListOfMolecules[i]->ScanForPeriodicCorrection(out);
 | 
|---|
 | 220 | 
 | 
|---|
 | 221 |     // output xyz file
 | 
|---|
 | 222 |     FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
 | 
|---|
 | 223 |     sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
 | 
|---|
 | 224 |     outputFragment.open(FragmentName, ios::out);
 | 
|---|
| [1e9384] | 225 |     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
 | 
|---|
| [db942e] | 226 |     if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
 | 
|---|
 | 227 |       *out << " done." << endl;
 | 
|---|
 | 228 |     else
 | 
|---|
 | 229 |       *out << " failed." << endl;
 | 
|---|
 | 230 |     result = result && intermediateResult;
 | 
|---|
 | 231 |     outputFragment.close();
 | 
|---|
 | 232 |     outputFragment.clear();
 | 
|---|
 | 233 | 
 | 
|---|
 | 234 |     *out << Verbose(2) << "Contained atoms: ";
 | 
|---|
 | 235 |     Walker = ListOfMolecules[i]->start;
 | 
|---|
 | 236 |     while (Walker->next != ListOfMolecules[i]->end) {
 | 
|---|
 | 237 |       Walker = Walker->next;
 | 
|---|
 | 238 |       *out << Walker->Name << " ";
 | 
|---|
| [14de469] | 239 |     }
 | 
|---|
| [db942e] | 240 |     *out << endl;
 | 
|---|
 | 241 |     // prepare output of config file
 | 
|---|
 | 242 |     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
 | 
|---|
 | 243 |     outputFragment.open(FragmentName, ios::out);
 | 
|---|
 | 244 |     strcpy(PathBackup, configuration->configpath);
 | 
|---|
 | 245 |     sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
 | 
|---|
 | 246 |     
 | 
|---|
 | 247 |     // center on edge
 | 
|---|
 | 248 |     ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
 | 
|---|
 | 249 |     ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
 | 
|---|
 | 250 |     int j = -1;
 | 
|---|
| [7f3b9d] | 251 |     for (int k=0;k<NDIM;k++) {
 | 
|---|
| [db942e] | 252 |       j += k+1;
 | 
|---|
 | 253 |       BoxDimension.x[k] = 5.;
 | 
|---|
 | 254 |       ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
 | 
|---|
| [14de469] | 255 |     }
 | 
|---|
| [db942e] | 256 |     ListOfMolecules[i]->Translate(&BoxDimension);
 | 
|---|
 | 257 | 
 | 
|---|
 | 258 |     // also calculate necessary orbitals
 | 
|---|
 | 259 |     ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
 | 
|---|
 | 260 |     ListOfMolecules[i]->CalculateOrbitals(*configuration);
 | 
|---|
 | 261 |     
 | 
|---|
 | 262 |     // change path in config
 | 
|---|
 | 263 |     configuration->SetDefaultPath(FragmentName);
 | 
|---|
| [1e9384] | 264 |     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
 | 
|---|
| [db942e] | 265 |     if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
 | 
|---|
 | 266 |       *out << " done." << endl;
 | 
|---|
 | 267 |     else
 | 
|---|
 | 268 |       *out << " failed." << endl;
 | 
|---|
 | 269 |     // restore old config
 | 
|---|
 | 270 |     configuration->SetDefaultPath(PathBackup);
 | 
|---|
 | 271 |     
 | 
|---|
 | 272 |     result = result && intermediateResult;
 | 
|---|
 | 273 |     outputFragment.close();
 | 
|---|
 | 274 |     outputFragment.clear();
 | 
|---|
| [7f3b9d] | 275 |     delete(FragmentNumber);
 | 
|---|
 | 276 |     //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
 | 
|---|
| [14de469] | 277 |   }
 | 
|---|
| [2459b1] | 278 |   cout << " done." << endl;
 | 
|---|
| [14de469] | 279 |   
 | 
|---|
 | 280 |   // printing final number
 | 
|---|
| [db942e] | 281 |   *out << "Final number of fragments: " << FragmentCounter << "." << endl;
 | 
|---|
| [14de469] | 282 |       
 | 
|---|
 | 283 |   return result;
 | 
|---|
 | 284 | };
 | 
|---|
 | 285 | 
 | 
|---|
 | 286 | /******************************************* Class MoleculeLeafClass ************************************************/
 | 
|---|
 | 287 | 
 | 
|---|
 | 288 | /** Constructor for MoleculeLeafClass root leaf.
 | 
|---|
 | 289 |  * \param *Up Leaf on upper level
 | 
|---|
 | 290 |  * \param *PreviousLeaf NULL - We are the first leaf on this level, otherwise points to previous in list
 | 
|---|
 | 291 |  */
 | 
|---|
 | 292 | //MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *Up = NULL, MoleculeLeafClass *Previous = NULL)
 | 
|---|
 | 293 | MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
 | 
|---|
 | 294 | {
 | 
|---|
 | 295 | //  if (Up != NULL)
 | 
|---|
 | 296 | //    if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
 | 
|---|
 | 297 | //      Up->DownLeaf = this;
 | 
|---|
 | 298 | //  UpLeaf = Up;
 | 
|---|
 | 299 | //  DownLeaf = NULL;
 | 
|---|
 | 300 |   Leaf = NULL;
 | 
|---|
 | 301 |   previous = PreviousLeaf;
 | 
|---|
 | 302 |   if (previous != NULL) {
 | 
|---|
 | 303 |     MoleculeLeafClass *Walker = previous->next;
 | 
|---|
 | 304 |     previous->next = this;
 | 
|---|
 | 305 |     next = Walker;
 | 
|---|
 | 306 |   } else {
 | 
|---|
 | 307 |     next = NULL;
 | 
|---|
 | 308 |   }
 | 
|---|
 | 309 | };
 | 
|---|
 | 310 | 
 | 
|---|
 | 311 | /** Destructor for MoleculeLeafClass.
 | 
|---|
 | 312 |  */
 | 
|---|
 | 313 | MoleculeLeafClass::~MoleculeLeafClass()
 | 
|---|
 | 314 | {
 | 
|---|
 | 315 | //  if (DownLeaf != NULL) {// drop leaves further down
 | 
|---|
 | 316 | //    MoleculeLeafClass *Walker = DownLeaf;
 | 
|---|
 | 317 | //    MoleculeLeafClass *Next;
 | 
|---|
 | 318 | //    do {
 | 
|---|
 | 319 | //      Next = Walker->NextLeaf;
 | 
|---|
 | 320 | //      delete(Walker);
 | 
|---|
 | 321 | //      Walker = Next;
 | 
|---|
 | 322 | //    } while (Walker != NULL);
 | 
|---|
 | 323 | //    // Last Walker sets DownLeaf automatically to NULL
 | 
|---|
 | 324 | //  }
 | 
|---|
 | 325 |   // remove the leaf itself
 | 
|---|
 | 326 |   if (Leaf != NULL) {
 | 
|---|
 | 327 |     delete(Leaf);
 | 
|---|
 | 328 |     Leaf = NULL;
 | 
|---|
 | 329 |   }
 | 
|---|
 | 330 |   // remove this Leaf from level list
 | 
|---|
 | 331 |   if (previous != NULL)   
 | 
|---|
 | 332 |     previous->next = next;
 | 
|---|
 | 333 | //  } else { // we are first in list (connects to UpLeaf->DownLeaf)
 | 
|---|
 | 334 | //    if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
 | 
|---|
 | 335 | //      NextLeaf->UpLeaf = UpLeaf;  // either null as we are top level or the upleaf of the first node
 | 
|---|
 | 336 | //    if (UpLeaf != NULL)
 | 
|---|
 | 337 | //      UpLeaf->DownLeaf = NextLeaf;  // either null as we are only leaf or NextLeaf if we are just the first
 | 
|---|
 | 338 | //  }
 | 
|---|
 | 339 | //  UpLeaf = NULL;
 | 
|---|
 | 340 |   if (next != NULL) // are we last in list
 | 
|---|
 | 341 |     next->previous = previous;
 | 
|---|
 | 342 |   next = NULL;
 | 
|---|
 | 343 |   previous = NULL;
 | 
|---|
 | 344 | };
 | 
|---|
 | 345 | 
 | 
|---|
 | 346 | /** Adds \a molecule leaf to the tree.
 | 
|---|
 | 347 |  * \param *ptr ptr to molecule to be added
 | 
|---|
 | 348 |  * \param *Previous previous MoleculeLeafClass referencing level and which on the level
 | 
|---|
 | 349 |  * \return true - success, false - something went wrong
 | 
|---|
 | 350 |  */
 | 
|---|
 | 351 | bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
 | 
|---|
 | 352 | {
 | 
|---|
 | 353 |   return false;
 | 
|---|
 | 354 | };
 | 
|---|
| [c82f3d] | 355 | 
 | 
|---|
 | 356 | /** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
 | 
|---|
 | 357 |  * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
 | 
|---|
 | 358 |  * \param *out output stream for debugging
 | 
|---|
 | 359 |  * \param *reference reference molecule with the bond structure to be copied
 | 
|---|
 | 360 |  * \param &FragmentCounter Counter needed to address \a **ListOfLocalAtoms
 | 
|---|
 | 361 |  * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
 | 
|---|
 | 362 |  * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 | 
|---|
 | 363 |  * \return true - success, false - faoilure
 | 
|---|
 | 364 |  */
 | 
|---|
 | 365 | bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
 | 
|---|
 | 366 | {
 | 
|---|
 | 367 |   atom *Walker = NULL, *OtherWalker = NULL;
 | 
|---|
 | 368 |   bond *Binder = NULL;
 | 
|---|
 | 369 |   bool status = true;
 | 
|---|
 | 370 |   int AtomNo;
 | 
|---|
 | 371 | 
 | 
|---|
| [d2a294] | 372 |   // fill ListOfLocalAtoms if NULL was given
 | 
|---|
 | 373 |   if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
 | 
|---|
 | 374 |     *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
 | 
|---|
 | 375 |     return false;
 | 
|---|
| [c82f3d] | 376 |   }
 | 
|---|
 | 377 |   
 | 
|---|
 | 378 |   if (status) {
 | 
|---|
 | 379 |     *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
 | 
|---|
 | 380 |     Walker = Leaf->start;
 | 
|---|
 | 381 |     while (Walker->next != Leaf->end) {
 | 
|---|
 | 382 |       Walker = Walker->next;
 | 
|---|
 | 383 |       AtomNo = Walker->father->nr;  // global id of the current walker
 | 
|---|
 | 384 |       for(int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
 | 
|---|
 | 385 |         Binder = reference->ListOfBondsPerAtom[AtomNo][i];
 | 
|---|
 | 386 |         OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->father)->nr];    // local copy of current bond partner of walker
 | 
|---|
 | 387 |         if (OtherWalker != NULL) {
 | 
|---|
 | 388 |           if (OtherWalker->nr > Walker->nr)
 | 
|---|
 | 389 |           Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
 | 
|---|
 | 390 |         } else {
 | 
|---|
 | 391 |           *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->father)->nr << "] is NULL!" << endl;
 | 
|---|
 | 392 |           status = false;
 | 
|---|
 | 393 |         }
 | 
|---|
 | 394 |       }
 | 
|---|
 | 395 |     }
 | 
|---|
 | 396 |     Leaf->CreateListOfBondsPerAtom(out);
 | 
|---|
 | 397 |     FragmentCounter++;
 | 
|---|
 | 398 |     if (next != NULL)
 | 
|---|
 | 399 |       status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
 | 
|---|
 | 400 |   }
 | 
|---|
 | 401 |   
 | 
|---|
 | 402 |   if (FreeList) {
 | 
|---|
 | 403 |     // free the index lookup list
 | 
|---|
| [386aa2] | 404 |     Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
 | 
|---|
 | 405 |     if (ListOfLocalAtoms[FragmentCounter] == NULL)
 | 
|---|
| [c82f3d] | 406 |       Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
 | 
|---|
 | 407 |   }
 | 
|---|
| [386aa2] | 408 |   FragmentCounter--;
 | 
|---|
| [c82f3d] | 409 |   return status;
 | 
|---|
 | 410 | };
 | 
|---|
 | 411 | 
 | 
|---|
| [386aa2] | 412 | /** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
 | 
|---|
 | 413 |  * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
 | 
|---|
 | 414 |  * \param *out output stream for debugging
 | 
|---|
 | 415 |  * \param *&RootStack stack to be filled
 | 
|---|
| [5de3c9] | 416 |  * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site
 | 
|---|
| [386aa2] | 417 |  * \param &FragmentCounter counts through the fragments in this MoleculeLeafClass
 | 
|---|
 | 418 |  * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
 | 
|---|
 | 419 |  */
 | 
|---|
| [de293ac] | 420 | bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
 | 
|---|
| [386aa2] | 421 | {
 | 
|---|
| [5de3c9] | 422 |   atom *Walker = NULL, *Father = NULL;
 | 
|---|
| [386aa2] | 423 | 
 | 
|---|
 | 424 |   if (RootStack != NULL) {
 | 
|---|
| [de293ac] | 425 |     // find first root candidates 
 | 
|---|
 | 426 |     if (&(RootStack[FragmentCounter]) != NULL) {
 | 
|---|
 | 427 |       RootStack[FragmentCounter].clear();  
 | 
|---|
 | 428 |       Walker = Leaf->start;
 | 
|---|
 | 429 |       while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
 | 
|---|
 | 430 |         Walker = Walker->next;
 | 
|---|
 | 431 |         Father = Walker->GetTrueFather();
 | 
|---|
 | 432 |         if (AtomMask[Father->nr]) // apply mask
 | 
|---|
 | 433 |     #ifdef ADDHYDROGEN
 | 
|---|
 | 434 |           if (Walker->type->Z != 1) // skip hydrogen
 | 
|---|
 | 435 |     #endif
 | 
|---|
 | 436 |             RootStack[FragmentCounter].push_front(Walker->nr);
 | 
|---|
| [d2a294] | 437 |       }
 | 
|---|
| [de293ac] | 438 |       if (next != NULL)
 | 
|---|
 | 439 |         next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
 | 
|---|
 | 440 |     }  else {
 | 
|---|
 | 441 |       *out << Verbose(1) << "Rootstack[" << FragmentCounter  << "] is NULL." << endl;
 | 
|---|
 | 442 |       return false;
 | 
|---|
| [386aa2] | 443 |     }
 | 
|---|
| [d2a294] | 444 |     FragmentCounter--;
 | 
|---|
| [386aa2] | 445 |     return true;
 | 
|---|
| [d2a294] | 446 |   } else {
 | 
|---|
 | 447 |     *out << Verbose(1) << "Rootstack is NULL." << endl;
 | 
|---|
| [386aa2] | 448 |     return false;
 | 
|---|
| [d2a294] | 449 |   }
 | 
|---|
 | 450 | };
 | 
|---|
 | 451 | 
 | 
|---|
 | 452 | /** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
 | 
|---|
 | 453 |  * \param *out output stream fro debugging
 | 
|---|
 | 454 |  * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 | 
|---|
 | 455 |  * \param &FragmentCounter counts the fragments as we move along the list
 | 
|---|
 | 456 |  * \param GlobalAtomCount number of atoms in the complete molecule
 | 
|---|
 | 457 |  * \param &FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 | 
|---|
 | 458 |  * \return true - succes, false - failure
 | 
|---|
 | 459 |  */
 | 
|---|
 | 460 | bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, int &FragmentCounter, int GlobalAtomCount, bool &FreeList)
 | 
|---|
 | 461 | {
 | 
|---|
 | 462 |   bool status = true;
 | 
|---|
 | 463 |   
 | 
|---|
 | 464 |   int Counter = Count();
 | 
|---|
 | 465 |   if (ListOfLocalAtoms == NULL) { // allocated initial pointer
 | 
|---|
 | 466 |     // allocate and set each field to NULL
 | 
|---|
 | 467 |     ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
 | 
|---|
 | 468 |     if (ListOfLocalAtoms != NULL) {
 | 
|---|
| [7f3b9d] | 469 |       for (int i=Counter;i--;)
 | 
|---|
| [d2a294] | 470 |         ListOfLocalAtoms[i] = NULL;
 | 
|---|
 | 471 |       FreeList = FreeList && true;
 | 
|---|
 | 472 |     } else
 | 
|---|
 | 473 |       status = false;
 | 
|---|
 | 474 |   }
 | 
|---|
 | 475 |   
 | 
|---|
 | 476 |   if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
 | 
|---|
 | 477 |     status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
 | 
|---|
 | 478 |     FreeList = FreeList && true;
 | 
|---|
 | 479 |   }
 | 
|---|
 | 480 |   
 | 
|---|
 | 481 |   return status;
 | 
|---|
 | 482 | };
 | 
|---|
 | 483 | 
 | 
|---|
 | 484 | /** The indices per keyset are compared to the respective father's Atom::nr in each subgraph and thus put into \a **&FragmentList.
 | 
|---|
 | 485 |  * \param *out output stream fro debugging
 | 
|---|
 | 486 |  * \param *reference reference molecule with the bond structure to be copied
 | 
|---|
 | 487 |  * \param *KeySetList list with all keysets
 | 
|---|
 | 488 |  * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 | 
|---|
 | 489 |  * \param **&FragmentList list to be allocated and returned
 | 
|---|
 | 490 |  * \param &FragmentCounter counts the fragments as we move along the list
 | 
|---|
 | 491 |  * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 | 
|---|
 | 492 |  * \retuen true - success, false - failure
 | 
|---|
 | 493 |  */
 | 
|---|
 | 494 | bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
 | 
|---|
 | 495 | {
 | 
|---|
 | 496 |   bool status = true;
 | 
|---|
 | 497 |   int KeySetCounter = 0;
 | 
|---|
 | 498 |   
 | 
|---|
 | 499 |   // fill ListOfLocalAtoms if NULL was given
 | 
|---|
 | 500 |   if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
 | 
|---|
 | 501 |     *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
 | 
|---|
 | 502 |     return false;
 | 
|---|
 | 503 |   }
 | 
|---|
 | 504 | 
 | 
|---|
 | 505 |   // allocate fragment list
 | 
|---|
 | 506 |   if (FragmentList == NULL) {
 | 
|---|
 | 507 |     KeySetCounter = Count();
 | 
|---|
| [bf46da] | 508 |     FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
 | 
|---|
| [7f3b9d] | 509 |     for(int i=KeySetCounter;i--;)
 | 
|---|
| [d2a294] | 510 |       FragmentList[i] = NULL;
 | 
|---|
 | 511 |     KeySetCounter = 0;
 | 
|---|
 | 512 |   }
 | 
|---|
 | 513 |   
 | 
|---|
 | 514 |   if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
 | 
|---|
 | 515 |     // assign scanned keysets
 | 
|---|
 | 516 |     if (FragmentList[FragmentCounter] == NULL)
 | 
|---|
 | 517 |       FragmentList[FragmentCounter] = new Graph;
 | 
|---|
 | 518 |     KeySet *TempSet = new KeySet;
 | 
|---|
 | 519 |     for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
 | 
|---|
 | 520 |       if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr]->nr != -1) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
 | 
|---|
 | 521 |         // translate keyset to local numbers
 | 
|---|
 | 522 |         for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
 | 
|---|
 | 523 |           TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
 | 
|---|
 | 524 |         // insert into FragmentList
 | 
|---|
 | 525 |         FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
 | 
|---|
 | 526 |       }
 | 
|---|
 | 527 |       TempSet->clear();
 | 
|---|
 | 528 |     }
 | 
|---|
 | 529 |     delete(TempSet);
 | 
|---|
 | 530 |     if (KeySetCounter == 0) {// if there are no keysets, delete the list
 | 
|---|
 | 531 |       *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
 | 
|---|
 | 532 |       delete(FragmentList[FragmentCounter]);
 | 
|---|
 | 533 |     } else
 | 
|---|
 | 534 |       *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
 | 
|---|
 | 535 |     FragmentCounter++;
 | 
|---|
 | 536 |     if (next != NULL)
 | 
|---|
 | 537 |       next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
 | 
|---|
 | 538 |     FragmentCounter--;
 | 
|---|
 | 539 |   } else
 | 
|---|
 | 540 |     *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
 | 
|---|
 | 541 |   
 | 
|---|
 | 542 |   return status;
 | 
|---|
| [386aa2] | 543 | };
 | 
|---|
| [c82f3d] | 544 | 
 | 
|---|
| [87f6c9] | 545 | /** Translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
 | 
|---|
 | 546 |  * \param *out output stream for debugging
 | 
|---|
 | 547 |  * \param **FragmentList Graph with local numbers per fragment
 | 
|---|
 | 548 |  * \param &FragmentCounter counts the fragments as we move along the list
 | 
|---|
 | 549 |  * \param &TotalNumberOfKeySets global key set counter
 | 
|---|
 | 550 |  * \param &TotalGraph Graph to be filled with global numbers
 | 
|---|
 | 551 |  */ 
 | 
|---|
 | 552 | void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
 | 
|---|
 | 553 | {
 | 
|---|
 | 554 |   KeySet *TempSet = new KeySet;
 | 
|---|
| [de293ac] | 555 |   if (FragmentList[FragmentCounter] != NULL) {
 | 
|---|
 | 556 |     for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
 | 
|---|
 | 557 |       for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
 | 
|---|
 | 558 |         TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
 | 
|---|
 | 559 |       TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
 | 
|---|
 | 560 |       TempSet->clear();
 | 
|---|
 | 561 |     }
 | 
|---|
 | 562 |     delete(TempSet);
 | 
|---|
 | 563 |   } else {
 | 
|---|
 | 564 |     *out << Verbose(1) << "FragmentList is NULL." << endl;
 | 
|---|
| [87f6c9] | 565 |   }
 | 
|---|
 | 566 |   if (next != NULL)
 | 
|---|
 | 567 |     next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
 | 
|---|
 | 568 |   FragmentCounter--;
 | 
|---|
 | 569 | };
 | 
|---|
 | 570 | 
 | 
|---|
| [386aa2] | 571 | /** Simply counts the number of items in the list, from given MoleculeLeafClass.
 | 
|---|
 | 572 |  * \return number of items
 | 
|---|
 | 573 |  */
 | 
|---|
| [d2a294] | 574 | int MoleculeLeafClass::Count() const
 | 
|---|
| [386aa2] | 575 | {
 | 
|---|
 | 576 |   if (next != NULL)
 | 
|---|
 | 577 |     return next->Count()+1;
 | 
|---|
 | 578 |   else
 | 
|---|
 | 579 |     return 1; 
 | 
|---|
 | 580 | };
 | 
|---|