Changeset d067d45 for src/moleculelist.cpp
- Timestamp:
- Jul 23, 2009, 1:45:24 PM (15 years ago)
- 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:
- 51c910
- Parents:
- ce5ac3 (diff), 437922 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the(diff)
links above to see all the changes relative to each parent. - git-author:
- Frederik Heber <heber@…> (07/23/09 12:34:47)
- git-committer:
- Frederik Heber <heber@…> (07/23/09 13:45:24)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/moleculelist.cpp
rce5ac3 rd067d45 1 1 /** \file MoleculeListClass.cpp 2 * 2 * 3 3 * Function implementations for the class MoleculeListClass. 4 * 4 * 5 5 */ 6 6 … … 13 13 MoleculeListClass::MoleculeListClass() 14 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"); 24 for (int i=NumMolecules;i--;) 25 ListOfMolecules[i] = NULL; 26 NumberOfMolecules = NumMolecules; 27 NumberOfTopAtoms = NumAtoms; 28 }; 29 15 // empty lists 16 ListOfMolecules.clear(); 17 MaxIndex = 1; 18 }; 30 19 31 20 /** Destructor for MoleculeListClass. … … 34 23 { 35 24 cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl; 36 for (int i=NumberOfMolecules;i--;) { 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 } 25 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) { 26 cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl; 27 delete (*ListRunner); 42 28 } 43 29 cout << Verbose(4) << "Freeing ListOfMolecules." << endl; 44 Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules"); 30 ListOfMolecules.clear(); // empty list 31 }; 32 33 /** Insert a new molecule into the list and set its number. 34 * \param *mol molecule to add to list. 35 * \return true - add successful 36 */ 37 void MoleculeListClass::insert(molecule *mol) 38 { 39 mol->IndexNr = MaxIndex++; 40 ListOfMolecules.push_back(mol); 45 41 }; 46 42 … … 57 53 atom *aWalker = NULL; 58 54 atom *bWalker = NULL; 59 55 60 56 // sort each atom list and put the numbers into a list, then go through 61 57 //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl; 62 if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount) {58 if ((**(molecule **) a).AtomCount < (**(molecule **) b).AtomCount) { 63 59 return -1; 64 } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount) 65 return +1; 60 } else { 61 if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount) 62 return +1; 66 63 else { 67 Count = (**(molecule **) a).AtomCount;64 Count = (**(molecule **) a).AtomCount; 68 65 aList = new int[Count]; 69 66 bList = new int[Count]; 70 67 71 68 // fill the lists 72 aWalker = (**(molecule **) a).start;73 bWalker = (**(molecule **) b).start;69 aWalker = (**(molecule **) a).start; 70 bWalker = (**(molecule **) b).start; 74 71 Counter = 0; 75 72 aCounter = 0; 76 73 bCounter = 0; 77 while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **)b).end)) {74 while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) { 78 75 aWalker = aWalker->next; 79 76 bWalker = bWalker->next; … … 90 87 // check if AtomCount was for real 91 88 flag = 0; 92 if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **)b).end)) {89 if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) { 93 90 flag = -1; 94 91 } else { 95 if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **)b).end))92 if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end)) 96 93 flag = 1; 97 94 } 98 95 if (flag == 0) { 99 96 // sort the lists 100 gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);101 gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);102 // compare the lists 103 97 gsl_heapsort(aList, Count, sizeof(int), CompareDoubles); 98 gsl_heapsort(bList, Count, sizeof(int), CompareDoubles); 99 // compare the lists 100 104 101 flag = 0; 105 for (int i=0;i<Count;i++) {102 for (int i = 0; i < Count; i++) { 106 103 if (aList[i] < bList[i]) { 107 104 flag = -1; … … 114 111 } 115 112 } 116 delete[] (aList);117 delete[] (bList);113 delete[] (aList); 114 delete[] (bList); 118 115 return flag; 119 116 } 120 117 } 121 return -1; 118 return -1; 119 }; 120 121 /** Output of a list of all molecules. 122 * \param *out output stream 123 */ 124 void MoleculeListClass::Enumerate(ofstream *out) 125 { 126 element* Elemental = NULL; 127 atom *Walker = NULL; 128 int Counts[MAX_ELEMENTS]; 129 double size=0; 130 Vector Origin; 131 132 // header 133 *out << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl; 134 cout << Verbose(0) << "-----------------------------------------------" << endl; 135 if (ListOfMolecules.size() == 0) 136 *out << "\tNone" << endl; 137 else { 138 Origin.Zero(); 139 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) { 140 // reset element counts 141 for (int j = 0; j<MAX_ELEMENTS;j++) 142 Counts[j] = 0; 143 // count atoms per element and determine size of bounding sphere 144 size=0.; 145 Walker = (*ListRunner)->start; 146 while (Walker->next != (*ListRunner)->end) { 147 Walker = Walker->next; 148 Counts[Walker->type->Z]++; 149 if (Walker->x.DistanceSquared(&Origin) > size) 150 size = Walker->x.DistanceSquared(&Origin); 151 } 152 // output Index, Name, number of atoms, chemical formula 153 *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t"; 154 Elemental = (*ListRunner)->elemente->end; 155 while(Elemental->previous != (*ListRunner)->elemente->start) { 156 Elemental = Elemental->previous; 157 if (Counts[Elemental->Z] != 0) 158 *out << Elemental->symbol << Counts[Elemental->Z]; 159 } 160 // Center and size 161 *out << "\t" << (*ListRunner)->Center << "\t" << sqrt(size) << endl; 162 } 163 } 164 }; 165 166 /** Returns the molecule with the given index \a index. 167 * \param index index of the desired molecule 168 * \return pointer to molecule structure, NULL if not found 169 */ 170 molecule * MoleculeListClass::ReturnIndex(int index) 171 { 172 for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) 173 if ((*ListRunner)->IndexNr == index) 174 return (*ListRunner); 175 return NULL; 176 }; 177 178 /** Simple merge of two molecules into one. 179 * \param *mol destination molecule 180 * \param *srcmol source molecule 181 * \return true - merge successful, false - merge failed (probably due to non-existant indices 182 */ 183 bool MoleculeListClass::SimpleMerge(molecule *mol, molecule *srcmol) 184 { 185 if (srcmol == NULL) 186 return false; 187 188 // put all molecules of src into mol 189 atom *Walker = srcmol->start; 190 atom *NextAtom = Walker->next; 191 while (NextAtom != srcmol->end) { 192 Walker = NextAtom; 193 NextAtom = Walker->next; 194 srcmol->UnlinkAtom(Walker); 195 mol->AddAtom(Walker); 196 } 197 198 // remove src 199 ListOfMolecules.remove(srcmol); 200 delete(srcmol); 201 return true; 202 }; 203 204 /** Simple add of one molecules into another. 205 * \param *mol destination molecule 206 * \param *srcmol source molecule 207 * \return true - merge successful, false - merge failed (probably due to non-existant indices 208 */ 209 bool MoleculeListClass::SimpleAdd(molecule *mol, molecule *srcmol) 210 { 211 if (srcmol == NULL) 212 return false; 213 214 // put all molecules of src into mol 215 atom *Walker = srcmol->start; 216 atom *NextAtom = Walker->next; 217 while (NextAtom != srcmol->end) { 218 Walker = NextAtom; 219 NextAtom = Walker->next; 220 Walker = mol->AddCopyAtom(Walker); 221 Walker->father = Walker; 222 } 223 224 return true; 225 }; 226 227 /** Simple merge of a given set of molecules into one. 228 * \param *mol destination molecule 229 * \param *src index of set of source molecule 230 * \param N number of source molecules 231 * \return true - merge successful, false - some merges failed (probably due to non-existant indices) 232 */ 233 bool MoleculeListClass::SimpleMultiMerge(molecule *mol, int *src, int N) 234 { 235 bool status = true; 236 // check presence of all source molecules 237 for (int i=0;i<N;i++) { 238 molecule *srcmol = ReturnIndex(src[i]); 239 status = status && SimpleMerge(mol, srcmol); 240 } 241 return status; 242 }; 243 244 /** Simple add of a given set of molecules into one. 245 * \param *mol destination molecule 246 * \param *src index of set of source molecule 247 * \param N number of source molecules 248 * \return true - merge successful, false - some merges failed (probably due to non-existant indices) 249 */ 250 bool MoleculeListClass::SimpleMultiAdd(molecule *mol, int *src, int N) 251 { 252 bool status = true; 253 // check presence of all source molecules 254 for (int i=0;i<N;i++) { 255 molecule *srcmol = ReturnIndex(src[i]); 256 status = status && SimpleAdd(mol, srcmol); 257 } 258 return status; 259 }; 260 261 /** Scatter merge of a given set of molecules into one. 262 * Scatter merge distributes the molecules in such a manner that they don't overlap. 263 * \param *mol destination molecule 264 * \param *src index of set of source molecule 265 * \param N number of source molecules 266 * \return true - merge successful, false - merge failed (probably due to non-existant indices 267 * \TODO find scatter center for each src molecule 268 */ 269 bool MoleculeListClass::ScatterMerge(molecule *mol, int *src, int N) 270 { 271 // check presence of all source molecules 272 for (int i=0;i<N;i++) { 273 // get pointer to src molecule 274 molecule *srcmol = ReturnIndex(src[i]); 275 if (srcmol == NULL) 276 return false; 277 } 278 // adapt each Center 279 for (int i=0;i<N;i++) { 280 // get pointer to src molecule 281 molecule *srcmol = ReturnIndex(src[i]); 282 //srcmol->Center.Zero(); 283 srcmol->Translate(&srcmol->Center); 284 } 285 // perform a simple multi merge 286 SimpleMultiMerge(mol, src, N); 287 return true; 288 }; 289 290 /** Embedding merge of a given set of molecules into one. 291 * Embedding merge inserts one molecule into the other. 292 * \param *mol destination molecule 293 * \param *srcmol source molecule 294 * \return true - merge successful, false - merge failed (probably due to non-existant indices 295 * \TODO find embedding center 296 */ 297 bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol) 298 { 299 if (srcmol == NULL) 300 return false; 301 302 // calculate center for merge 303 srcmol->Center.CopyVector(mol->FindEmbeddingHole((ofstream *)&cout, srcmol)); 304 srcmol->Center.Zero(); 305 306 // perform simple merge 307 SimpleMerge(mol, srcmol); 308 return true; 122 309 }; 123 310 … … 127 314 void MoleculeListClass::Output(ofstream *out) 128 315 { 129 *out << Verbose(1) << "MoleculeList: ";130 for ( int i=0;i<NumberOfMolecules;i++)131 *out << ListOfMolecules[i]<< "\t";316 *out << Verbose(1) << "MoleculeList: "; 317 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) 318 *out << *ListRunner << "\t"; 132 319 *out << endl; 133 320 }; … … 145 332 atom *Runner = NULL; 146 333 double ***FitConstant = NULL, **correction = NULL; 147 int a, b;334 int a, b; 148 335 ofstream output; 149 336 ifstream input; … … 165 352 input.open(line.c_str()); 166 353 if (input == NULL) { 167 cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl; 354 cerr << endl << "Unable to open " << line << ", is the directory correct?" 355 << endl; 168 356 return false; 169 357 } 170 a =0;171 b =-1; // we overcount by one358 a = 0; 359 b = -1; // we overcount by one 172 360 while (!input.eof()) { 173 361 input.getline(ParsedLine, 1023); 174 362 zeile.str(ParsedLine); 175 int i =0;363 int i = 0; 176 364 while (!zeile.eof()) { 177 365 zeile >> distance; 178 366 i++; 179 367 } 180 if (i > a) 181 a = i; 368 if (i > a) 369 a = i; 182 370 b++; 183 371 } 184 372 cout << "I recognized " << a << " columns and " << b << " rows, "; 185 373 input.close(); 186 374 187 375 // 0b. allocate memory for constants 188 FitConstant = (double ***) Malloc(sizeof(double **) *3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");189 for (int k =0;k<3;k++) {190 FitConstant[k] = (double **) Malloc(sizeof(double *) *a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");191 for (int i =a;i--;) {192 FitConstant[k][i] = (double *) Malloc(sizeof(double) *b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");376 FitConstant = (double ***) Malloc(sizeof(double **) * 3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant"); 377 for (int k = 0; k < 3; k++) { 378 FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]"); 379 for (int i = a; i--;) { 380 FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]"); 193 381 } 194 382 } 195 383 // 0c. parse in constants 196 for (int i =0;i<3;i++) {384 for (int i = 0; i < 3; i++) { 197 385 line = path; 198 386 line.append("/"); 199 387 line += FRAGMENTPREFIX; 200 sprintf(ParsedLine, "%d", i +1);388 sprintf(ParsedLine, "%d", i + 1); 201 389 line += ParsedLine; 202 390 line += FITCONSTANTSUFFIX; … … 206 394 return false; 207 395 } 208 int k = 0, l;396 int k = 0, l; 209 397 while ((!input.eof()) && (k < b)) { 210 398 input.getline(ParsedLine, 1023); … … 223 411 input.close(); 224 412 } 225 for (int k=0;k<3;k++) {413 for (int k = 0; k < 3; k++) { 226 414 cout << "Constants " << k << ":" << endl; 227 for (int j =0;j<b;j++) {228 for (int i =0;i<a;i++) {415 for (int j = 0; j < b; j++) { 416 for (int i = 0; i < a; i++) { 229 417 cout << FitConstant[k][i][j] << "\t"; 230 418 } … … 233 421 cout << endl; 234 422 } 235 423 236 424 // 0d. allocate final correction matrix 237 correction = (double **) Malloc(sizeof(double *) *a, "MoleculeListClass::AddHydrogenCorrection: **correction");238 for (int i =a;i--;)239 correction[i] = (double *) Malloc(sizeof(double) *b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");240 425 correction = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **correction"); 426 for (int i = a; i--;) 427 correction[i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *correction[]"); 428 241 429 // 1a. go through every molecule in the list 242 for (int i=NumberOfMolecules;i--;) {430 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) { 243 431 // 1b. zero final correction matrix 244 for (int k =a;k--;)245 for (int j =b;j--;)432 for (int k = a; k--;) 433 for (int j = b; j--;) 246 434 correction[k][j] = 0.; 247 435 // 2. take every hydrogen that is a saturated one 248 Walker = ListOfMolecules[i]->start;249 while (Walker->next != ListOfMolecules[i]->end) {436 Walker = (*ListRunner)->start; 437 while (Walker->next != (*ListRunner)->end) { 250 438 Walker = Walker->next; 251 //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl; 252 if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen 253 Runner = ListOfMolecules[i]->start; 254 while (Runner->next != ListOfMolecules[i]->end) { 439 //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl; 440 if ((Walker->type->Z == 1) && ((Walker->father == NULL) 441 || (Walker->father->type->Z != 1))) { // if it's a hydrogen 442 Runner = (*ListRunner)->start; 443 while (Runner->next != (*ListRunner)->end) { 255 444 Runner = Runner->next; 256 //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << * ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;445 //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl; 257 446 // 3. take every other hydrogen that is the not the first and not bound to same bonding partner 258 if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ( ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {// (hydrogens have only one bonding partner!)447 if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!) 259 448 // 4. evaluate the morse potential for each matrix component and add up 260 distance = sqrt(Runner->x.Distance(&Walker->x));261 //cout << "Fragment " << i<< ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;262 for (int k=0;k<a;k++) {263 for (int j =0;j<b;j++) {264 switch (k) {265 case 1:266 case 7:267 case 11:268 tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);269 break;270 default:271 tmp = FitConstant[0][k][j] * pow( distance,FitConstant[1][k][j]) + FitConstant[2][k][j];449 distance = Runner->x.Distance(&Walker->x); 450 //cout << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl; 451 for (int k = 0; k < a; k++) { 452 for (int j = 0; j < b; j++) { 453 switch (k) { 454 case 1: 455 case 7: 456 case 11: 457 tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2); 458 break; 459 default: 460 tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j]; 272 461 }; 273 correction[k][j] -= tmp; 462 correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction) 274 463 //cout << tmp << "\t"; 275 464 } … … 285 474 line.append("/"); 286 475 line += FRAGMENTPREFIX; 287 FragmentNumber = FixedDigitNumber( NumberOfMolecules, i);476 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr); 288 477 line += FragmentNumber; 289 delete (FragmentNumber);478 delete (FragmentNumber); 290 479 line += HCORRECTIONSUFFIX; 291 480 output.open(line.c_str()); 292 481 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl; 293 for (int j =0;j<b;j++) {294 for (int i=0;i<a;i++)482 for (int j = 0; j < b; j++) { 483 for (int i = 0; i < a; i++) 295 484 output << correction[i][j] << "\t"; 296 485 output << endl; … … 303 492 output.open(line.c_str()); 304 493 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl; 305 for (int j =0;j<b;j++) {306 for (int i=0;i<a;i++)494 for (int j = 0; j < b; j++) { 495 for (int i = 0; i < a; i++) 307 496 output << 0 << "\t"; 308 497 output << endl; … … 310 499 output.close(); 311 500 // 6. free memory of parsed matrices 312 FitConstant = (double ***) Malloc(sizeof(double **) *a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");313 for (int k =0;k<3;k++) {314 FitConstant[k] = (double **) Malloc(sizeof(double *) *a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");315 for (int i =a;i--;) {316 FitConstant[k][i] = (double *) Malloc(sizeof(double) *b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");501 FitConstant = (double ***) Malloc(sizeof(double **) * a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant"); 502 for (int k = 0; k < 3; k++) { 503 FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]"); 504 for (int i = a; i--;) { 505 FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]"); 317 506 } 318 507 } … … 327 516 * \return true - file written successfully, false - writing failed 328 517 */ 329 bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex) 518 bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, 519 int *SortIndex) 330 520 { 331 521 bool status = true; … … 342 532 //cout << Verbose(1) << "Final AtomicForcesList: "; 343 533 //output << prefix << "Forces" << endl; 344 for(int j=0;j<NumberOfMolecules;j++) { 345 //if (TEList[j] != 0) { 346 runner = ListOfMolecules[j]->elemente->start; 347 while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element 534 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) { 535 runner = (*ListRunner)->elemente->start; 536 while (runner->next != (*ListRunner)->elemente->end) { // go through every element 348 537 runner = runner->next; 349 if ( ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms350 Walker = ListOfMolecules[j]->start;351 while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element538 if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms 539 Walker = (*ListRunner)->start; 540 while (Walker->next != (*ListRunner)->end) { // go through every atom of this element 352 541 Walker = Walker->next; 353 542 if (Walker->type->Z == runner->Z) { 354 543 if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea 355 544 //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it 356 ForcesFile << 357 } else // otherwise a -1 to indicate an added saturation hydrogen358 ForcesFile << "-1\t";359 }545 ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t"; 546 } else 547 // otherwise a -1 to indicate an added saturation hydrogen 548 ForcesFile << "-1\t"; 360 549 } 361 550 } 551 } 362 552 } 363 553 ForcesFile << endl; … … 370 560 } 371 561 ForcesFile.close(); 372 562 373 563 return status; 374 564 }; … … 382 572 * \return true - success (each file was written), false - something went wrong. 383 573 */ 384 bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, con st char *fragmentprefix, config *configuration, int *SortIndex, bool DoPeriodic, bool DoCentering)574 bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex) 385 575 { 386 576 ofstream outputFragment; … … 395 585 int FragmentCounter = 0; 396 586 ofstream output; 397 string basis("3-21G"); 398 587 399 588 // store the fragments as config and as xyz 400 for (int i=0;i<NumberOfMolecules;i++) {589 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) { 401 590 // save default path as it is changed for each fragment 402 591 path = configuration->GetDefaultPath(); … … 407 596 408 597 // correct periodic 409 if (DoPeriodic) 410 ListOfMolecules[i]->ScanForPeriodicCorrection(out); 598 (*ListRunner)->ScanForPeriodicCorrection(out); 411 599 412 600 // output xyz file 413 FragmentNumber = FixedDigitNumber( NumberOfMolecules, FragmentCounter++);414 sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, fragmentprefix, FragmentNumber);601 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++); 602 sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber); 415 603 outputFragment.open(FragmentName, ios::out); 416 *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";417 if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))604 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ..."; 605 if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment))) 418 606 *out << " done." << endl; 419 607 else … … 425 613 // list atoms in fragment for debugging 426 614 *out << Verbose(2) << "Contained atoms: "; 427 Walker = ListOfMolecules[i]->start;428 while (Walker->next != ListOfMolecules[i]->end) {615 Walker = (*ListRunner)->start; 616 while (Walker->next != (*ListRunner)->end) { 429 617 Walker = Walker->next; 430 618 *out << Walker->Name << " "; 431 619 } 432 620 *out << endl; 433 621 434 622 // center on edge 435 if (DoCentering) { 436 ListOfMolecules[i]->CenterEdge(out, &BoxDimension); 437 ListOfMolecules[i]->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary 438 int j = -1; 439 for (int k=0;k<NDIM;k++) { 440 j += k+1; 441 BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem); 442 ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.; 443 } 444 ListOfMolecules[i]->Translate(&BoxDimension); 445 } 623 (*ListRunner)->CenterEdge(out, &BoxDimension); 624 (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary 625 int j = -1; 626 for (int k = 0; k < NDIM; k++) { 627 j += k + 1; 628 BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem); 629 (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.; 630 } 631 (*ListRunner)->Translate(&BoxDimension); 446 632 447 633 // also calculate necessary orbitals 448 ListOfMolecules[i]->CountElements(); // this is a bugfix, atoms should should actually be added correctly to this fragment449 ListOfMolecules[i]->CalculateOrbitals(*configuration);450 634 (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment 635 (*ListRunner)->CalculateOrbitals(*configuration); 636 451 637 // change path in config 452 638 //strcpy(PathBackup, configuration->configpath); 453 sprintf(FragmentName, "%s/%s%s/", PathBackup, fragmentprefix, FragmentNumber);639 sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber); 454 640 configuration->SetDefaultPath(FragmentName); 455 641 456 642 // and save as config 457 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, fragmentprefix, FragmentNumber);458 *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";459 if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))643 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber); 644 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ..."; 645 if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner)))) 460 646 *out << " done." << endl; 461 647 else … … 466 652 configuration->SetDefaultPath(PathBackup); 467 653 468 469 654 // and save as mpqc input file 470 sprintf(FragmentName, "%s/%s%s. in", configuration->configpath, fragmentprefix, FragmentNumber);471 *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";472 if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))655 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber); 656 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ..."; 657 if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner)))) 473 658 *out << " done." << endl; 474 659 else 475 660 *out << " failed." << endl; 476 661 477 662 result = result && intermediateResult; 478 663 //outputFragment.close(); 479 664 //outputFragment.clear(); 480 delete (FragmentNumber);665 delete (FragmentNumber); 481 666 //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber"); 482 667 } 483 668 cout << " done." << endl; 484 669 485 670 // printing final number 486 671 *out << "Final number of fragments: " << FragmentCounter << "." << endl; 487 672 488 673 return result; 489 674 }; 675 676 /** Counts the number of molecules with the molecule::ActiveFlag set. 677 * \return number of molecules with ActiveFlag set to true. 678 */ 679 int MoleculeListClass::NumberOfActiveMolecules() 680 { 681 int count = 0; 682 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) 683 count += ((*ListRunner)->ActiveFlag ? 1 : 0); 684 return count; 685 }; 686 490 687 491 688 /******************************************* Class MoleculeLeafClass ************************************************/ … … 498 695 MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL) 499 696 { 500 // if (Up != NULL)501 // if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?502 // Up->DownLeaf = this;503 // UpLeaf = Up;504 // DownLeaf = NULL;697 // if (Up != NULL) 698 // if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf? 699 // Up->DownLeaf = this; 700 // UpLeaf = Up; 701 // DownLeaf = NULL; 505 702 Leaf = NULL; 506 703 previous = PreviousLeaf; … … 518 715 MoleculeLeafClass::~MoleculeLeafClass() 519 716 { 520 // if (DownLeaf != NULL) {// drop leaves further down521 // MoleculeLeafClass *Walker = DownLeaf;522 // MoleculeLeafClass *Next;523 // do {524 // Next = Walker->NextLeaf;525 // delete(Walker);526 // Walker = Next;527 // } while (Walker != NULL);528 // // Last Walker sets DownLeaf automatically to NULL529 // }717 // if (DownLeaf != NULL) {// drop leaves further down 718 // MoleculeLeafClass *Walker = DownLeaf; 719 // MoleculeLeafClass *Next; 720 // do { 721 // Next = Walker->NextLeaf; 722 // delete(Walker); 723 // Walker = Next; 724 // } while (Walker != NULL); 725 // // Last Walker sets DownLeaf automatically to NULL 726 // } 530 727 // remove the leaf itself 531 728 if (Leaf != NULL) { 532 delete (Leaf);729 delete (Leaf); 533 730 Leaf = NULL; 534 731 } 535 732 // remove this Leaf from level list 536 if (previous != NULL) 733 if (previous != NULL) 537 734 previous->next = next; 538 // } else { // we are first in list (connects to UpLeaf->DownLeaf)539 // if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))540 // NextLeaf->UpLeaf = UpLeaf; // either null as we are top level or the upleaf of the first node541 // if (UpLeaf != NULL)542 // UpLeaf->DownLeaf = NextLeaf; // either null as we are only leaf or NextLeaf if we are just the first543 // }544 // UpLeaf = NULL;735 // } else { // we are first in list (connects to UpLeaf->DownLeaf) 736 // if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL)) 737 // NextLeaf->UpLeaf = UpLeaf; // either null as we are top level or the upleaf of the first node 738 // if (UpLeaf != NULL) 739 // UpLeaf->DownLeaf = NextLeaf; // either null as we are only leaf or NextLeaf if we are just the first 740 // } 741 // UpLeaf = NULL; 545 742 if (next != NULL) // are we last in list 546 743 next->previous = previous; … … 581 778 return false; 582 779 } 583 780 584 781 if (status) { 585 *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl; 782 *out << Verbose(1) << "Creating adjacency list for subgraph " << this 783 << "." << endl; 586 784 Walker = Leaf->start; 587 785 while (Walker->next != Leaf->end) { 588 786 Walker = Walker->next; 589 AtomNo = Walker->GetTrueFather()->nr; 590 for (int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all787 AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker 788 for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all 591 789 Binder = reference->ListOfBondsPerAtom[AtomNo][i]; 592 OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; 790 OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker 593 791 if (OtherWalker != NULL) { 594 792 if (OtherWalker->nr > Walker->nr) 595 Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);793 Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree); 596 794 } else { 597 795 *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl; … … 606 804 FragmentCounter--; 607 805 } 608 806 609 807 if ((FreeList) && (ListOfLocalAtoms != NULL)) { 610 808 // free the index lookup list 611 Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");809 Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]"); 612 810 if (FragmentCounter == 0) // first fragments frees the initial pointer to list 613 Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");811 Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms"); 614 812 } 615 813 FragmentCounter--; … … 626 824 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update 627 825 */ 628 bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter) 826 bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, 827 KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter) 629 828 { 630 829 atom *Walker = NULL, *Father = NULL; 631 830 632 831 if (RootStack != NULL) { 633 // find first root candidates 832 // find first root candidates 634 833 if (&(RootStack[FragmentCounter]) != NULL) { 635 RootStack[FragmentCounter].clear(); 834 RootStack[FragmentCounter].clear(); 636 835 Walker = Leaf->start; 637 836 while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms … … 639 838 Father = Walker->GetTrueFather(); 640 839 if (AtomMask[Father->nr]) // apply mask 641 840 #ifdef ADDHYDROGEN 642 841 if (Walker->type->Z != 1) // skip hydrogen 643 644 842 #endif 843 RootStack[FragmentCounter].push_front(Walker->nr); 645 844 } 646 845 if (next != NULL) 647 846 next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter); 648 } 649 *out << Verbose(1) << "Rootstack[" << FragmentCounter 847 } else { 848 *out << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl; 650 849 return false; 651 850 } … … 669 868 { 670 869 bool status = true; 671 870 672 871 int Counter = Count(); 673 872 if (ListOfLocalAtoms == NULL) { // allocated initial pointer 674 873 // allocate and set each field to NULL 675 ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **) *Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");874 ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **) * Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms"); 676 875 if (ListOfLocalAtoms != NULL) { 677 for (int i =Counter;i--;)876 for (int i = Counter; i--;) 678 877 ListOfLocalAtoms[i] = NULL; 679 878 FreeList = FreeList && true; … … 681 880 status = false; 682 881 } 683 882 684 883 if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph 685 884 status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount); 686 885 FreeList = FreeList && true; 687 886 } 688 887 689 888 return status; 690 889 }; … … 700 899 * \retuen true - success, false - failure 701 900 */ 702 bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList) 901 bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, 902 molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, 903 Graph **&FragmentList, int &FragmentCounter, bool FreeList) 703 904 { 704 905 bool status = true; 705 906 int KeySetCounter = 0; 706 907 707 908 *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl; 708 909 // fill ListOfLocalAtoms if NULL was given … … 715 916 if (FragmentList == NULL) { 716 917 KeySetCounter = Count(); 717 FragmentList = (Graph **) Malloc(sizeof(Graph *) *KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");718 for (int i=KeySetCounter;i--;)918 FragmentList = (Graph **) Malloc(sizeof(Graph *) * KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList"); 919 for (int i = KeySetCounter; i--;) 719 920 FragmentList[i] = NULL; 720 921 KeySetCounter = 0; 721 922 } 722 923 723 924 if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all 724 925 // assign scanned keysets … … 726 927 FragmentList[FragmentCounter] = new Graph; 727 928 KeySet *TempSet = new KeySet; 728 for (Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!729 if ( 929 for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers! 930 if (ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set 730 931 // translate keyset to local numbers 731 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)932 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++) 732 933 TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr); 733 934 // insert into FragmentList 734 FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));935 FragmentList[FragmentCounter]->insert(GraphPair(*TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second))); 735 936 } 736 937 TempSet->clear(); 737 938 } 738 delete (TempSet);939 delete (TempSet); 739 940 if (KeySetCounter == 0) {// if there are no keysets, delete the list 740 941 *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl; 741 delete (FragmentList[FragmentCounter]);942 delete (FragmentList[FragmentCounter]); 742 943 } else 743 944 *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl; … … 748 949 } else 749 950 *out << Verbose(1) << "KeySetList is NULL or empty." << endl; 750 951 751 952 if ((FreeList) && (ListOfLocalAtoms != NULL)) { 752 953 // free the index lookup list 753 Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");954 Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]"); 754 955 if (FragmentCounter == 0) // first fragments frees the initial pointer to list 755 Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");956 Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms"); 756 957 } 757 958 *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl; … … 765 966 * \param &TotalNumberOfKeySets global key set counter 766 967 * \param &TotalGraph Graph to be filled with global numbers 767 */ 768 void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph) 968 */ 969 void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, 970 Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, 971 Graph &TotalGraph) 769 972 { 770 973 *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl; 771 974 KeySet *TempSet = new KeySet; 772 975 if (FragmentList[FragmentCounter] != NULL) { 773 for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {774 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)976 for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) { 977 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++) 775 978 TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr); 776 TotalGraph.insert(GraphPair(*TempSet, pair<int, double>(TotalNumberOfKeySets++, (*runner).second.second)));979 TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second))); 777 980 TempSet->clear(); 778 981 } 779 delete (TempSet);982 delete (TempSet); 780 983 } else { 781 984 *out << Verbose(1) << "FragmentList is NULL." << endl; … … 793 996 { 794 997 if (next != NULL) 795 return next->Count() +1;998 return next->Count() + 1; 796 999 else 797 return 1; 798 }; 1000 return 1; 1001 }; 1002
Note:
See TracChangeset
for help on using the changeset viewer.