source: src/moleculelist.cpp@ 7ea9e6

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
Last change on this file since 7ea9e6 was c9bce3e, checked in by Frederik Heber <heber@…>, 16 years ago

config::Load() refactored: Dissection into connected subgraphs -> MoleculeListClass::DissectMoleculeIntoConnectedSubgraphs()

  • Property mode set to 100755
File size: 42.8 KB
Line 
1/** \file MoleculeListClass.cpp
2 *
3 * Function implementations for the class MoleculeListClass.
4 *
5 */
6
7#include "atom.hpp"
8#include "bond.hpp"
9#include "boundary.hpp"
10#include "config.hpp"
11#include "element.hpp"
12#include "helpers.hpp"
13#include "linkedcell.hpp"
14#include "lists.hpp"
15#include "molecule.hpp"
16#include "memoryallocator.hpp"
17#include "periodentafel.hpp"
18
19/*********************************** Functions for class MoleculeListClass *************************/
20
21/** Constructor for MoleculeListClass.
22 */
23MoleculeListClass::MoleculeListClass()
24{
25 // empty lists
26 ListOfMolecules.clear();
27 MaxIndex = 1;
28};
29
30/** Destructor for MoleculeListClass.
31 */
32MoleculeListClass::~MoleculeListClass()
33{
34 cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
35 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
36 cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl;
37 delete (*ListRunner);
38 }
39 cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
40 ListOfMolecules.clear(); // empty list
41};
42
43/** Insert a new molecule into the list and set its number.
44 * \param *mol molecule to add to list.
45 * \return true - add successful
46 */
47void MoleculeListClass::insert(molecule *mol)
48{
49 mol->IndexNr = MaxIndex++;
50 ListOfMolecules.push_back(mol);
51};
52
53/** Compare whether two molecules are equal.
54 * \param *a molecule one
55 * \param *n molecule two
56 * \return lexical value (-1, 0, +1)
57 */
58int MolCompare(const void *a, const void *b)
59{
60 int *aList = NULL, *bList = NULL;
61 int Count, Counter, aCounter, bCounter;
62 int flag;
63 atom *aWalker = NULL;
64 atom *bWalker = NULL;
65
66 // sort each atom list and put the numbers into a list, then go through
67 //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
68 if ((**(molecule **) a).AtomCount < (**(molecule **) b).AtomCount) {
69 return -1;
70 } else {
71 if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount)
72 return +1;
73 else {
74 Count = (**(molecule **) a).AtomCount;
75 aList = new int[Count];
76 bList = new int[Count];
77
78 // fill the lists
79 aWalker = (**(molecule **) a).start;
80 bWalker = (**(molecule **) b).start;
81 Counter = 0;
82 aCounter = 0;
83 bCounter = 0;
84 while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
85 aWalker = aWalker->next;
86 bWalker = bWalker->next;
87 if (aWalker->GetTrueFather() == NULL)
88 aList[Counter] = Count + (aCounter++);
89 else
90 aList[Counter] = aWalker->GetTrueFather()->nr;
91 if (bWalker->GetTrueFather() == NULL)
92 bList[Counter] = Count + (bCounter++);
93 else
94 bList[Counter] = bWalker->GetTrueFather()->nr;
95 Counter++;
96 }
97 // check if AtomCount was for real
98 flag = 0;
99 if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
100 flag = -1;
101 } else {
102 if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end))
103 flag = 1;
104 }
105 if (flag == 0) {
106 // sort the lists
107 gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
108 gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
109 // compare the lists
110
111 flag = 0;
112 for (int i = 0; i < Count; i++) {
113 if (aList[i] < bList[i]) {
114 flag = -1;
115 } else {
116 if (aList[i] > bList[i])
117 flag = 1;
118 }
119 if (flag != 0)
120 break;
121 }
122 }
123 delete[] (aList);
124 delete[] (bList);
125 return flag;
126 }
127 }
128 return -1;
129};
130
131/** Output of a list of all molecules.
132 * \param *out output stream
133 */
134void MoleculeListClass::Enumerate(ofstream *out)
135{
136 element* Elemental = NULL;
137 atom *Walker = NULL;
138 int Counts[MAX_ELEMENTS];
139 double size=0;
140 Vector Origin;
141
142 // header
143 *out << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
144 cout << Verbose(0) << "-----------------------------------------------" << endl;
145 if (ListOfMolecules.size() == 0)
146 *out << "\tNone" << endl;
147 else {
148 Origin.Zero();
149 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
150 // reset element counts
151 for (int j = 0; j<MAX_ELEMENTS;j++)
152 Counts[j] = 0;
153 // count atoms per element and determine size of bounding sphere
154 size=0.;
155 Walker = (*ListRunner)->start;
156 while (Walker->next != (*ListRunner)->end) {
157 Walker = Walker->next;
158 Counts[Walker->type->Z]++;
159 if (Walker->x.DistanceSquared(&Origin) > size)
160 size = Walker->x.DistanceSquared(&Origin);
161 }
162 // output Index, Name, number of atoms, chemical formula
163 *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
164 Elemental = (*ListRunner)->elemente->end;
165 while(Elemental->previous != (*ListRunner)->elemente->start) {
166 Elemental = Elemental->previous;
167 if (Counts[Elemental->Z] != 0)
168 *out << Elemental->symbol << Counts[Elemental->Z];
169 }
170 // Center and size
171 *out << "\t" << (*ListRunner)->Center << "\t" << sqrt(size) << endl;
172 }
173 }
174};
175
176/** Returns the molecule with the given index \a index.
177 * \param index index of the desired molecule
178 * \return pointer to molecule structure, NULL if not found
179 */
180molecule * MoleculeListClass::ReturnIndex(int index)
181{
182 for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
183 if ((*ListRunner)->IndexNr == index)
184 return (*ListRunner);
185 return NULL;
186};
187
188/** Simple merge of two molecules into one.
189 * \param *mol destination molecule
190 * \param *srcmol source molecule
191 * \return true - merge successful, false - merge failed (probably due to non-existant indices
192 */
193bool MoleculeListClass::SimpleMerge(molecule *mol, molecule *srcmol)
194{
195 if (srcmol == NULL)
196 return false;
197
198 // put all molecules of src into mol
199 atom *Walker = srcmol->start;
200 atom *NextAtom = Walker->next;
201 while (NextAtom != srcmol->end) {
202 Walker = NextAtom;
203 NextAtom = Walker->next;
204 srcmol->UnlinkAtom(Walker);
205 mol->AddAtom(Walker);
206 }
207
208 // remove src
209 ListOfMolecules.remove(srcmol);
210 delete(srcmol);
211 return true;
212};
213
214/** Simple add of one molecules into another.
215 * \param *mol destination molecule
216 * \param *srcmol source molecule
217 * \return true - merge successful, false - merge failed (probably due to non-existant indices
218 */
219bool MoleculeListClass::SimpleAdd(molecule *mol, molecule *srcmol)
220{
221 if (srcmol == NULL)
222 return false;
223
224 // put all molecules of src into mol
225 atom *Walker = srcmol->start;
226 atom *NextAtom = Walker->next;
227 while (NextAtom != srcmol->end) {
228 Walker = NextAtom;
229 NextAtom = Walker->next;
230 Walker = mol->AddCopyAtom(Walker);
231 Walker->father = Walker;
232 }
233
234 return true;
235};
236
237/** Simple merge of a given set of molecules into one.
238 * \param *mol destination molecule
239 * \param *src index of set of source molecule
240 * \param N number of source molecules
241 * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
242 */
243bool MoleculeListClass::SimpleMultiMerge(molecule *mol, int *src, int N)
244{
245 bool status = true;
246 // check presence of all source molecules
247 for (int i=0;i<N;i++) {
248 molecule *srcmol = ReturnIndex(src[i]);
249 status = status && SimpleMerge(mol, srcmol);
250 }
251 return status;
252};
253
254/** Simple add of a given set of molecules into one.
255 * \param *mol destination molecule
256 * \param *src index of set of source molecule
257 * \param N number of source molecules
258 * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
259 */
260bool MoleculeListClass::SimpleMultiAdd(molecule *mol, int *src, int N)
261{
262 bool status = true;
263 // check presence of all source molecules
264 for (int i=0;i<N;i++) {
265 molecule *srcmol = ReturnIndex(src[i]);
266 status = status && SimpleAdd(mol, srcmol);
267 }
268 return status;
269};
270
271/** Scatter merge of a given set of molecules into one.
272 * Scatter merge distributes the molecules in such a manner that they don't overlap.
273 * \param *mol destination molecule
274 * \param *src index of set of source molecule
275 * \param N number of source molecules
276 * \return true - merge successful, false - merge failed (probably due to non-existant indices
277 * \TODO find scatter center for each src molecule
278 */
279bool MoleculeListClass::ScatterMerge(molecule *mol, int *src, int N)
280{
281 // check presence of all source molecules
282 for (int i=0;i<N;i++) {
283 // get pointer to src molecule
284 molecule *srcmol = ReturnIndex(src[i]);
285 if (srcmol == NULL)
286 return false;
287 }
288 // adapt each Center
289 for (int i=0;i<N;i++) {
290 // get pointer to src molecule
291 molecule *srcmol = ReturnIndex(src[i]);
292 //srcmol->Center.Zero();
293 srcmol->Translate(&srcmol->Center);
294 }
295 // perform a simple multi merge
296 SimpleMultiMerge(mol, src, N);
297 return true;
298};
299
300/** Embedding merge of a given set of molecules into one.
301 * Embedding merge inserts one molecule into the other.
302 * \param *mol destination molecule (fixed one)
303 * \param *srcmol source molecule (variable one, where atoms are taken from)
304 * \return true - merge successful, false - merge failed (probably due to non-existant indices)
305 * \TODO linked cell dimensions for boundary points has to be as big as inner diameter!
306 */
307bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
308{
309 LinkedCell *LCList = NULL;
310 Tesselation *TesselStruct = NULL;
311 if ((srcmol == NULL) || (mol == NULL)) {
312 cout << Verbose(1) << "ERROR: Either fixed or variable molecule is given as NULL." << endl;
313 return false;
314 }
315
316 // calculate envelope for *mol
317 LCList = new LinkedCell(mol, 8.);
318 FindNonConvexBorder((ofstream *)&cout, mol, TesselStruct, (const LinkedCell *&)LCList, 4., NULL);
319 if (TesselStruct == NULL) {
320 cout << Verbose(1) << "ERROR: Could not tesselate the fixed molecule." << endl;
321 return false;
322 }
323 delete(LCList);
324 LCList = new LinkedCell(TesselStruct, 8.); // re-create with boundary points only!
325
326 // prepare index list for bonds
327 srcmol->CountAtoms((ofstream *)&cout);
328 atom ** CopyAtoms = new atom*[srcmol->AtomCount];
329 for(int i=0;i<srcmol->AtomCount;i++)
330 CopyAtoms[i] = NULL;
331
332 // for each of the source atoms check whether we are in- or outside and add copy atom
333 atom *Walker = srcmol->start;
334 int nr=0;
335 while (Walker->next != srcmol->end) {
336 Walker = Walker->next;
337 cout << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl;
338 if (!TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
339 CopyAtoms[Walker->nr] = new atom(Walker);
340 mol->AddAtom(CopyAtoms[Walker->nr]);
341 nr++;
342 } else {
343 // do nothing
344 }
345 }
346 cout << Verbose(1) << nr << " of " << srcmol->AtomCount << " atoms have been merged.";
347
348 // go through all bonds and add as well
349 bond *Binder = srcmol->first;
350 while(Binder->next != srcmol->last) {
351 Binder = Binder->next;
352 cout << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
353 mol->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
354 }
355 delete(LCList);
356 return true;
357};
358
359/** Simple output of the pointers in ListOfMolecules.
360 * \param *out output stream
361 */
362void MoleculeListClass::Output(ofstream *out)
363{
364 *out << Verbose(1) << "MoleculeList: ";
365 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
366 *out << *ListRunner << "\t";
367 *out << endl;
368};
369
370/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
371 * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
372 * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
373 * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
374 * \param *out output stream for debugging
375 * \param *path path to file
376 */
377bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
378{
379 atom *Walker = NULL;
380 atom *Runner = NULL;
381 bond *Binder = NULL;
382 double ***FitConstant = NULL, **correction = NULL;
383 int a, b;
384 ofstream output;
385 ifstream input;
386 string line;
387 stringstream zeile;
388 double distance;
389 char ParsedLine[1023];
390 double tmp;
391 char *FragmentNumber = NULL;
392
393 cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
394 // 0. parse in fit constant files that should have the same dimension as the final energy files
395 // 0a. find dimension of matrices with constants
396 line = path;
397 line.append("/");
398 line += FRAGMENTPREFIX;
399 line += "1";
400 line += FITCONSTANTSUFFIX;
401 input.open(line.c_str());
402 if (input == NULL) {
403 cerr << endl << "Unable to open " << line << ", is the directory correct?"
404 << endl;
405 return false;
406 }
407 a = 0;
408 b = -1; // we overcount by one
409 while (!input.eof()) {
410 input.getline(ParsedLine, 1023);
411 zeile.str(ParsedLine);
412 int i = 0;
413 while (!zeile.eof()) {
414 zeile >> distance;
415 i++;
416 }
417 if (i > a)
418 a = i;
419 b++;
420 }
421 cout << "I recognized " << a << " columns and " << b << " rows, ";
422 input.close();
423
424 // 0b. allocate memory for constants
425 FitConstant = Calloc<double**>(3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
426 for (int k = 0; k < 3; k++) {
427 FitConstant[k] = Calloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
428 for (int i = a; i--;) {
429 FitConstant[k][i] = Calloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
430 }
431 }
432 // 0c. parse in constants
433 for (int i = 0; i < 3; i++) {
434 line = path;
435 line.append("/");
436 line += FRAGMENTPREFIX;
437 sprintf(ParsedLine, "%d", i + 1);
438 line += ParsedLine;
439 line += FITCONSTANTSUFFIX;
440 input.open(line.c_str());
441 if (input == NULL) {
442 cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
443 return false;
444 }
445 int k = 0, l;
446 while ((!input.eof()) && (k < b)) {
447 input.getline(ParsedLine, 1023);
448 //cout << "Current Line: " << ParsedLine << endl;
449 zeile.str(ParsedLine);
450 zeile.clear();
451 l = 0;
452 while ((!zeile.eof()) && (l < a)) {
453 zeile >> FitConstant[i][l][k];
454 //cout << FitConstant[i][l][k] << "\t";
455 l++;
456 }
457 //cout << endl;
458 k++;
459 }
460 input.close();
461 }
462 for (int k = 0; k < 3; k++) {
463 cout << "Constants " << k << ":" << endl;
464 for (int j = 0; j < b; j++) {
465 for (int i = 0; i < a; i++) {
466 cout << FitConstant[k][i][j] << "\t";
467 }
468 cout << endl;
469 }
470 cout << endl;
471 }
472
473 // 0d. allocate final correction matrix
474 correction = Calloc<double*>(a, "MoleculeListClass::AddHydrogenCorrection: **correction");
475 for (int i = a; i--;)
476 correction[i] = Calloc<double>(b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
477
478 // 1a. go through every molecule in the list
479 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
480 // 1b. zero final correction matrix
481 for (int k = a; k--;)
482 for (int j = b; j--;)
483 correction[k][j] = 0.;
484 // 2. take every hydrogen that is a saturated one
485 Walker = (*ListRunner)->start;
486 while (Walker->next != (*ListRunner)->end) {
487 Walker = Walker->next;
488 //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
489 if ((Walker->type->Z == 1) && ((Walker->father == NULL)
490 || (Walker->father->type->Z != 1))) { // if it's a hydrogen
491 Runner = (*ListRunner)->start;
492 while (Runner->next != (*ListRunner)->end) {
493 Runner = Runner->next;
494 //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
495 // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
496 Binder = *(Runner->ListOfBonds.begin());
497 if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (Binder->GetOtherAtom(Runner) != Binder->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
498 // 4. evaluate the morse potential for each matrix component and add up
499 distance = Runner->x.Distance(&Walker->x);
500 //cout << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
501 for (int k = 0; k < a; k++) {
502 for (int j = 0; j < b; j++) {
503 switch (k) {
504 case 1:
505 case 7:
506 case 11:
507 tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
508 break;
509 default:
510 tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
511 };
512 correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
513 //cout << tmp << "\t";
514 }
515 //cout << endl;
516 }
517 //cout << endl;
518 }
519 }
520 }
521 }
522 // 5. write final matrix to file
523 line = path;
524 line.append("/");
525 line += FRAGMENTPREFIX;
526 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
527 line += FragmentNumber;
528 delete (FragmentNumber);
529 line += HCORRECTIONSUFFIX;
530 output.open(line.c_str());
531 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
532 for (int j = 0; j < b; j++) {
533 for (int i = 0; i < a; i++)
534 output << correction[i][j] << "\t";
535 output << endl;
536 }
537 output.close();
538 }
539 line = path;
540 line.append("/");
541 line += HCORRECTIONSUFFIX;
542 output.open(line.c_str());
543 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
544 for (int j = 0; j < b; j++) {
545 for (int i = 0; i < a; i++)
546 output << 0 << "\t";
547 output << endl;
548 }
549 output.close();
550 // 6. free memory of parsed matrices
551 for (int k = 0; k < 3; k++) {
552 for (int i = a; i--;) {
553 Free(&FitConstant[k][i]);
554 }
555 Free(&FitConstant[k]);
556 }
557 Free(&FitConstant);
558 cout << "done." << endl;
559 return true;
560};
561
562/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
563 * \param *out output stream for debugging
564 * \param *path path to file
565 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
566 * \return true - file written successfully, false - writing failed
567 */
568bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path,
569 int *SortIndex)
570{
571 bool status = true;
572 ofstream ForcesFile;
573 stringstream line;
574 atom *Walker = NULL;
575 element *runner = NULL;
576
577 // open file for the force factors
578 *out << Verbose(1) << "Saving force factors ... ";
579 line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
580 ForcesFile.open(line.str().c_str(), ios::out);
581 if (ForcesFile != NULL) {
582 //cout << Verbose(1) << "Final AtomicForcesList: ";
583 //output << prefix << "Forces" << endl;
584 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
585 runner = (*ListRunner)->elemente->start;
586 while (runner->next != (*ListRunner)->elemente->end) { // go through every element
587 runner = runner->next;
588 if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms
589 Walker = (*ListRunner)->start;
590 while (Walker->next != (*ListRunner)->end) { // go through every atom of this element
591 Walker = Walker->next;
592 if (Walker->type->Z == runner->Z) {
593 if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
594 //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
595 ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t";
596 } else
597 // otherwise a -1 to indicate an added saturation hydrogen
598 ForcesFile << "-1\t";
599 }
600 }
601 }
602 }
603 ForcesFile << endl;
604 }
605 ForcesFile.close();
606 *out << Verbose(1) << "done." << endl;
607 } else {
608 status = false;
609 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
610 }
611 ForcesFile.close();
612
613 return status;
614};
615
616/** Writes a config file for each molecule in the given \a **FragmentList.
617 * \param *out output stream for debugging
618 * \param *configuration standard configuration to attach atoms in fragment molecule to.
619 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
620 * \param DoPeriodic true - call ScanForPeriodicCorrection, false - don't
621 * \param DoCentering true - call molecule::CenterEdge(), false - don't
622 * \return true - success (each file was written), false - something went wrong.
623 */
624bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
625{
626 ofstream outputFragment;
627 char FragmentName[MAXSTRINGSIZE];
628 char PathBackup[MAXSTRINGSIZE];
629 bool result = true;
630 bool intermediateResult = true;
631 atom *Walker = NULL;
632 Vector BoxDimension;
633 char *FragmentNumber = NULL;
634 char *path = NULL;
635 int FragmentCounter = 0;
636 ofstream output;
637
638 // store the fragments as config and as xyz
639 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
640 // save default path as it is changed for each fragment
641 path = configuration->GetDefaultPath();
642 if (path != NULL)
643 strcpy(PathBackup, path);
644 else
645 cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
646
647 // correct periodic
648 (*ListRunner)->ScanForPeriodicCorrection(out);
649
650 // output xyz file
651 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
652 sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
653 outputFragment.open(FragmentName, ios::out);
654 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ...";
655 if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
656 *out << " done." << endl;
657 else
658 *out << " failed." << endl;
659 result = result && intermediateResult;
660 outputFragment.close();
661 outputFragment.clear();
662
663 // list atoms in fragment for debugging
664 *out << Verbose(2) << "Contained atoms: ";
665 Walker = (*ListRunner)->start;
666 while (Walker->next != (*ListRunner)->end) {
667 Walker = Walker->next;
668 *out << Walker->Name << " ";
669 }
670 *out << endl;
671
672 // center on edge
673 (*ListRunner)->CenterEdge(out, &BoxDimension);
674 (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
675 int j = -1;
676 for (int k = 0; k < NDIM; k++) {
677 j += k + 1;
678 BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
679 (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
680 }
681 (*ListRunner)->Translate(&BoxDimension);
682
683 // also calculate necessary orbitals
684 (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
685 (*ListRunner)->CalculateOrbitals(*configuration);
686
687 // change path in config
688 //strcpy(PathBackup, configuration->configpath);
689 sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
690 configuration->SetDefaultPath(FragmentName);
691
692 // and save as config
693 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
694 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...";
695 if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner))))
696 *out << " done." << endl;
697 else
698 *out << " failed." << endl;
699 result = result && intermediateResult;
700
701 // restore old config
702 configuration->SetDefaultPath(PathBackup);
703
704 // and save as mpqc input file
705 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
706 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...";
707 if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner))))
708 *out << " done." << endl;
709 else
710 *out << " failed." << endl;
711
712 result = result && intermediateResult;
713 //outputFragment.close();
714 //outputFragment.clear();
715 Free(&FragmentNumber);
716 }
717 cout << " done." << endl;
718
719 // printing final number
720 *out << "Final number of fragments: " << FragmentCounter << "." << endl;
721
722 return result;
723};
724
725/** Counts the number of molecules with the molecule::ActiveFlag set.
726 * \return number of molecules with ActiveFlag set to true.
727 */
728int MoleculeListClass::NumberOfActiveMolecules()
729{
730 int count = 0;
731 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
732 count += ((*ListRunner)->ActiveFlag ? 1 : 0);
733 return count;
734};
735
736/** Dissects given \a *mol into connected subgraphs and inserts them as new molecules but with old atoms into \a this.
737 * \param *out output stream for debugging
738 * \param *mol molecule with atoms to dissect
739 * \param *configuration config with BondGraph
740 */
741void MoleculeListClass::DissectMoleculeIntoConnectedSubgraphs(ofstream * const out, molecule * const mol, config * const configuration)
742{
743 // 1. dissect the molecule into connected subgraphs
744 configuration ->BG->ConstructBondGraph(out, mol);
745
746 // 2. scan for connected subgraphs
747 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
748 class StackClass<bond *> *BackEdgeStack = NULL;
749 Subgraphs = mol->DepthFirstSearchAnalysis(out, BackEdgeStack);
750 delete(BackEdgeStack);
751
752 // 3. dissect (the following construct is needed to have the atoms not in the order of the DFS, but in
753 // the original one as parsed in)
754 // TODO: Optimize this, when molecules just contain pointer list of global atoms!
755
756 // 4a. create array of molecules to fill
757 const int MolCount = Subgraphs->next->Count();
758 molecule **molecules = Malloc<molecule *>(MolCount, "config::Load() - **molecules");
759 for (int i=0;i<MolCount;i++) {
760 molecules[i] = (molecule*) new molecule(mol->elemente);
761 molecules[i]->ActiveFlag = true;
762 insert(molecules[i]);
763 }
764
765 // 4b. create and fill map of which atom is associated to which connected molecule (note, counting starts at 1)
766 int FragmentCounter = 0;
767 int *MolMap = Calloc<int>(mol->AtomCount, "config::Load() - *MolMap");
768 MoleculeLeafClass *MolecularWalker = Subgraphs;
769 atom *Walker = NULL;
770 while (MolecularWalker->next != NULL) {
771 MolecularWalker = MolecularWalker->next;
772 Walker = MolecularWalker->Leaf->start;
773 while (Walker->next != MolecularWalker->Leaf->end) {
774 Walker = Walker->next;
775 MolMap[Walker->GetTrueFather()->nr] = FragmentCounter+1;
776 }
777 FragmentCounter++;
778 }
779
780 // 4c. relocate atoms to new molecules and remove from Leafs
781 Walker = mol->start;
782 while (mol->start->next != mol->end) {
783 Walker = mol->start->next;
784 if ((Walker->nr <0) || (Walker->nr >= mol->AtomCount)) {
785 cerr << "Index of atom " << *Walker << " is invalid!" << endl;
786 performCriticalExit();
787 }
788 FragmentCounter = MolMap[Walker->nr];
789 if (FragmentCounter != 0) {
790 cout << Verbose(3) << "Re-linking " << *Walker << "..." << endl;
791 unlink(Walker);
792 molecules[FragmentCounter-1]->AddAtom(Walker); // counting starts at 1
793 } else {
794 cerr << "Atom " << *Walker << " not associated to molecule!" << endl;
795 performCriticalExit();
796 }
797 }
798 // 4d. we don't need to redo bonds, as they are connected subgraphs and still maintained their ListOfBonds
799 // 4e. free Leafs
800 MolecularWalker = Subgraphs;
801 while (MolecularWalker->next != NULL) {
802 MolecularWalker = MolecularWalker->next;
803 delete(MolecularWalker->previous);
804 }
805 delete(MolecularWalker);
806 Free(&MolMap);
807 Free(&molecules);
808 cout << Verbose(1) << "I scanned " << FragmentCounter << " molecules." << endl;
809};
810
811/******************************************* Class MoleculeLeafClass ************************************************/
812
813/** Constructor for MoleculeLeafClass root leaf.
814 * \param *Up Leaf on upper level
815 * \param *PreviousLeaf NULL - We are the first leaf on this level, otherwise points to previous in list
816 */
817//MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *Up = NULL, MoleculeLeafClass *Previous = NULL)
818MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
819{
820 // if (Up != NULL)
821 // if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
822 // Up->DownLeaf = this;
823 // UpLeaf = Up;
824 // DownLeaf = NULL;
825 Leaf = NULL;
826 previous = PreviousLeaf;
827 if (previous != NULL) {
828 MoleculeLeafClass *Walker = previous->next;
829 previous->next = this;
830 next = Walker;
831 } else {
832 next = NULL;
833 }
834};
835
836/** Destructor for MoleculeLeafClass.
837 */
838MoleculeLeafClass::~MoleculeLeafClass()
839{
840 // if (DownLeaf != NULL) {// drop leaves further down
841 // MoleculeLeafClass *Walker = DownLeaf;
842 // MoleculeLeafClass *Next;
843 // do {
844 // Next = Walker->NextLeaf;
845 // delete(Walker);
846 // Walker = Next;
847 // } while (Walker != NULL);
848 // // Last Walker sets DownLeaf automatically to NULL
849 // }
850 // remove the leaf itself
851 if (Leaf != NULL) {
852 delete (Leaf);
853 Leaf = NULL;
854 }
855 // remove this Leaf from level list
856 if (previous != NULL)
857 previous->next = next;
858 // } else { // we are first in list (connects to UpLeaf->DownLeaf)
859 // if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
860 // NextLeaf->UpLeaf = UpLeaf; // either null as we are top level or the upleaf of the first node
861 // if (UpLeaf != NULL)
862 // UpLeaf->DownLeaf = NextLeaf; // either null as we are only leaf or NextLeaf if we are just the first
863 // }
864 // UpLeaf = NULL;
865 if (next != NULL) // are we last in list
866 next->previous = previous;
867 next = NULL;
868 previous = NULL;
869};
870
871/** Adds \a molecule leaf to the tree.
872 * \param *ptr ptr to molecule to be added
873 * \param *Previous previous MoleculeLeafClass referencing level and which on the level
874 * \return true - success, false - something went wrong
875 */
876bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
877{
878 return false;
879};
880
881/** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
882 * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
883 * \param *out output stream for debugging
884 * \param *reference reference molecule with the bond structure to be copied
885 * \param &FragmentCounter Counter needed to address \a **ListOfLocalAtoms
886 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
887 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
888 * \return true - success, false - faoilure
889 */
890bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, const molecule * const reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
891{
892 atom *Walker = NULL;
893 atom *OtherWalker = NULL;
894 atom *Father = NULL;
895 bool status = true;
896 int AtomNo;
897
898 *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
899 // fill ListOfLocalAtoms if NULL was given
900 if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
901 *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
902 return false;
903 }
904
905 if (status) {
906 *out << Verbose(1) << "Creating adjacency list for subgraph " << Leaf << "." << endl;
907 // remove every bond from the list
908 bond *Binder = NULL;
909 while (Leaf->last->previous != Leaf->first) {
910 Binder = Leaf->last->previous;
911 Binder->leftatom->UnregisterBond(Binder);
912 Binder->rightatom->UnregisterBond(Binder);
913 removewithoutcheck(Binder);
914 }
915
916 Walker = Leaf->start;
917 while (Walker->next != Leaf->end) {
918 Walker = Walker->next;
919 Father = Walker->GetTrueFather();
920 AtomNo = Father->nr; // global id of the current walker
921 for (BondList::const_iterator Runner = Father->ListOfBonds.begin(); Runner != Father->ListOfBonds.end(); (++Runner)) {
922 OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
923 if (OtherWalker != NULL) {
924 if (OtherWalker->nr > Walker->nr)
925 Leaf->AddBond(Walker, OtherWalker, (*Runner)->BondDegree);
926 } else {
927 *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
928 status = false;
929 }
930 }
931 }
932 }
933
934 if ((FreeList) && (ListOfLocalAtoms != NULL)) {
935 // free the index lookup list
936 Free(&ListOfLocalAtoms[FragmentCounter]);
937 if (FragmentCounter == 0) // first fragments frees the initial pointer to list
938 Free(&ListOfLocalAtoms);
939 }
940 *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
941 return status;
942};
943
944/** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
945 * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
946 * \param *out output stream for debugging
947 * \param *&RootStack stack to be filled
948 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site
949 * \param &FragmentCounter counts through the fragments in this MoleculeLeafClass
950 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
951 */
952bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out,
953 KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
954{
955 atom *Walker = NULL, *Father = NULL;
956
957 if (RootStack != NULL) {
958 // find first root candidates
959 if (&(RootStack[FragmentCounter]) != NULL) {
960 RootStack[FragmentCounter].clear();
961 Walker = Leaf->start;
962 while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
963 Walker = Walker->next;
964 Father = Walker->GetTrueFather();
965 if (AtomMask[Father->nr]) // apply mask
966#ifdef ADDHYDROGEN
967 if (Walker->type->Z != 1) // skip hydrogen
968#endif
969 RootStack[FragmentCounter].push_front(Walker->nr);
970 }
971 if (next != NULL)
972 next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
973 } else {
974 *out << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl;
975 return false;
976 }
977 FragmentCounter--;
978 return true;
979 } else {
980 *out << Verbose(1) << "Rootstack is NULL." << endl;
981 return false;
982 }
983};
984
985/** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
986 * \param *out output stream from debugging
987 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
988 * \param FragmentCounter counts the fragments as we move along the list
989 * \param GlobalAtomCount number of atoms in the complete molecule
990 * \param &FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
991 * \return true - success, false - failure
992 */
993bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList)
994{
995 bool status = true;
996
997 if (ListOfLocalAtoms == NULL) { // allocated initial pointer
998 // allocate and set each field to NULL
999 const int Counter = Count();
1000 ListOfLocalAtoms = Calloc<atom**>(Counter, "MoleculeLeafClass::FillListOfLocalAtoms - ***ListOfLocalAtoms");
1001 if (ListOfLocalAtoms == NULL) {
1002 FreeList = FreeList && false;
1003 status = false;
1004 }
1005 }
1006
1007 if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
1008 status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
1009 FreeList = FreeList && true;
1010 }
1011
1012 return status;
1013};
1014
1015/** The indices per keyset are compared to the respective father's Atom::nr in each subgraph and thus put into \a **&FragmentList.
1016 * \param *out output stream fro debugging
1017 * \param *reference reference molecule with the bond structure to be copied
1018 * \param *KeySetList list with all keysets
1019 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
1020 * \param **&FragmentList list to be allocated and returned
1021 * \param &FragmentCounter counts the fragments as we move along the list
1022 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
1023 * \retuen true - success, false - failure
1024 */
1025bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out,
1026 molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms,
1027 Graph **&FragmentList, int &FragmentCounter, bool FreeList)
1028{
1029 bool status = true;
1030 int KeySetCounter = 0;
1031
1032 *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
1033 // fill ListOfLocalAtoms if NULL was given
1034 if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
1035 *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
1036 return false;
1037 }
1038
1039 // allocate fragment list
1040 if (FragmentList == NULL) {
1041 KeySetCounter = Count();
1042 FragmentList = Calloc<Graph*>(KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
1043 KeySetCounter = 0;
1044 }
1045
1046 if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
1047 // assign scanned keysets
1048 if (FragmentList[FragmentCounter] == NULL)
1049 FragmentList[FragmentCounter] = new Graph;
1050 KeySet *TempSet = new KeySet;
1051 for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers!
1052 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
1053 // translate keyset to local numbers
1054 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
1055 TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
1056 // insert into FragmentList
1057 FragmentList[FragmentCounter]->insert(GraphPair(*TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
1058 }
1059 TempSet->clear();
1060 }
1061 delete (TempSet);
1062 if (KeySetCounter == 0) {// if there are no keysets, delete the list
1063 *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
1064 delete (FragmentList[FragmentCounter]);
1065 } else
1066 *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
1067 FragmentCounter++;
1068 if (next != NULL)
1069 next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
1070 FragmentCounter--;
1071 } else
1072 *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
1073
1074 if ((FreeList) && (ListOfLocalAtoms != NULL)) {
1075 // free the index lookup list
1076 Free(&ListOfLocalAtoms[FragmentCounter]);
1077 if (FragmentCounter == 0) // first fragments frees the initial pointer to list
1078 Free(&ListOfLocalAtoms);
1079 }
1080 *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
1081 return status;
1082};
1083
1084/** Translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
1085 * \param *out output stream for debugging
1086 * \param **FragmentList Graph with local numbers per fragment
1087 * \param &FragmentCounter counts the fragments as we move along the list
1088 * \param &TotalNumberOfKeySets global key set counter
1089 * \param &TotalGraph Graph to be filled with global numbers
1090 */
1091void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out,
1092 Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets,
1093 Graph &TotalGraph)
1094{
1095 *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
1096 KeySet *TempSet = new KeySet;
1097 if (FragmentList[FragmentCounter] != NULL) {
1098 for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
1099 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
1100 TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
1101 TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
1102 TempSet->clear();
1103 }
1104 delete (TempSet);
1105 } else {
1106 *out << Verbose(1) << "FragmentList is NULL." << endl;
1107 }
1108 if (next != NULL)
1109 next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
1110 FragmentCounter--;
1111 *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
1112};
1113
1114/** Simply counts the number of items in the list, from given MoleculeLeafClass.
1115 * \return number of items
1116 */
1117int MoleculeLeafClass::Count() const
1118{
1119 if (next != NULL)
1120 return next->Count() + 1;
1121 else
1122 return 1;
1123};
1124
Note: See TracBrowser for help on using the repository browser.