source: molecuilder/src/molecules.cpp@ e9c14d

Last change on this file since e9c14d was cc6dfb, checked in by Frederik Heber <heber@…>, 17 years ago

CyclicStructureAnalysis(): forget to check whether we actually have found rings or not (only then MinimumRingSize needs to be checked for non-set neighbours also)

  • Property mode set to 100644
File size: 185.3 KB
Line 
1/** \file molecules.cpp
2 *
3 * Functions for the class molecule.
4 *
5 */
6
7#include "molecules.hpp"
8
9/************************************* Other Functions *************************************/
10
11/** Determines sum of squared distances of \a X to all \a **vectors.
12 * \param *x reference vector
13 * \param *params
14 * \return sum of square distances
15 */
16double LSQ (const gsl_vector * x, void * params)
17{
18 double sum = 0.;
19 struct LSQ_params *par = (struct LSQ_params *)params;
20 vector **vectors = par->vectors;
21 int num = par->num;
22
23 for (int i=0;i<num;i++) {
24 for(int j=0;j<NDIM;j++)
25 sum += (gsl_vector_get(x,j) - (vectors[i])->x[j])*(gsl_vector_get(x,j) - (vectors[i])->x[j]);
26 }
27
28 return sum;
29};
30
31/************************************* Functions for class molecule *********************************/
32
33/** Constructor of class molecule.
34 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
35 */
36molecule::molecule(periodentafel *teil)
37{
38 // init atom chain list
39 start = new atom;
40 end = new atom;
41 start->father = NULL;
42 end->father = NULL;
43 link(start,end);
44 // init bond chain list
45 first = new bond(start, end, 1, -1);
46 last = new bond(start, end, 1, -1);
47 link(first,last);
48 // other stuff
49 last_atom = 0;
50 elemente = teil;
51 AtomCount = 0;
52 BondCount = 0;
53 NoNonBonds = 0;
54 NoNonHydrogen = 0;
55 NoCyclicBonds = 0;
56 ListOfBondsPerAtom = NULL;
57 NumberOfBondsPerAtom = NULL;
58 ElementCount = 0;
59 for(int i=0;i<MAX_ELEMENTS;i++)
60 ElementsInMolecule[i] = 0;
61 cell_size[0] = cell_size[2] = cell_size[5]= 20.;
62 cell_size[1] = cell_size[3] = cell_size[4]= 0.;
63};
64
65/** Destructor of class molecule.
66 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
67 */
68molecule::~molecule()
69{
70 if (ListOfBondsPerAtom != NULL)
71 for(int i=0;i<AtomCount;i++)
72 Free((void **)&ListOfBondsPerAtom[i], "molecule::~molecule: ListOfBondsPerAtom[i]");
73 Free((void **)&ListOfBondsPerAtom, "molecule::~molecule: ListOfBondsPerAtom");
74 Free((void **)&NumberOfBondsPerAtom, "molecule::~molecule: NumberOfBondsPerAtom");
75 CleanupMolecule();
76 delete(first);
77 delete(last);
78 delete(end);
79 delete(start);
80};
81
82/** Adds given atom \a *pointer from molecule list.
83 * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
84 * \param *pointer allocated and set atom
85 * \return true - succeeded, false - atom not found in list
86 */
87bool molecule::AddAtom(atom *pointer)
88{
89 if (pointer != NULL) {
90 pointer->sort = &pointer->nr;
91 pointer->nr = last_atom++; // increase number within molecule
92 AtomCount++;
93 if (pointer->type != NULL) {
94 if (ElementsInMolecule[pointer->type->Z] == 0)
95 ElementCount++;
96 ElementsInMolecule[pointer->type->Z]++; // increase number of elements
97 if (pointer->type->Z != 1)
98 NoNonHydrogen++;
99 if (pointer->Name == NULL) {
100 Free((void **)&pointer->Name, "molecule::AddAtom: *pointer->Name");
101 pointer->Name = (char *) Malloc(sizeof(char)*6, "molecule::AddAtom: *pointer->Name");
102 sprintf(pointer->Name, "%2s%02d", pointer->type->symbol, pointer->nr+1);
103 }
104 }
105 return add(pointer, end);
106 } else
107 return false;
108};
109
110/** Adds a copy of the given atom \a *pointer from molecule list.
111 * Increases molecule::last_atom and gives last number to added atom.
112 * \param *pointer allocated and set atom
113 * \return true - succeeded, false - atom not found in list
114 */
115atom * molecule::AddCopyAtom(atom *pointer)
116{
117 if (pointer != NULL) {
118 atom *walker = new atom();
119 walker->type = pointer->type; // copy element of atom
120 walker->x.CopyVector(&pointer->x); // copy coordination
121 walker->v.CopyVector(&pointer->v); // copy velocity
122 walker->FixedIon = pointer->FixedIon;
123 walker->sort = &walker->nr;
124 walker->nr = last_atom++; // increase number within molecule
125 walker->father = pointer; //->GetTrueFather();
126 walker->Name = (char *) Malloc(sizeof(char)*strlen(pointer->Name)+1, "molecule::AddCopyAtom: *Name");
127 strcpy (walker->Name, pointer->Name);
128 add(walker, end);
129 if ((pointer->type != NULL) && (pointer->type->Z != 1))
130 NoNonHydrogen++;
131 AtomCount++;
132 return walker;
133 } else
134 return NULL;
135};
136
137/** Adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
138 * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
139 * a different scheme when adding \a *replacement atom for the given one.
140 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
141 * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
142 * *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalVector().
143 * The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
144 * replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
145 * element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
146 * hydrogens forming this angle with *origin.
147 * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
148 * triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
149 * determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
150 * We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
151 * \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
152 * \f]
153 * vector::GetNormalVector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalVector creates
154 * the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
155 * The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
156 * the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
157 * \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
158 * \f]
159 * as the coordination of all three atoms in the coordinate system of these three vectors:
160 * \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
161 *
162 * \param *out output stream for debugging
163 * \param *Bond pointer to bond between \a *origin and \a *replacement
164 * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
165 * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
166 * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
167 * \param **BondList list of bonds \a *replacement has (necessary to determine plane for double and triple bonds)
168 * \param NumBond number of bonds in \a **BondList
169 * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
170 * \return number of atoms added, if < bond::BondDegree then something went wrong
171 * \todo double and triple bonds splitting (always use the tetraeder angle!)
172 */
173bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem)
174{
175 double bondlength; // bond length of the bond to be replaced/cut
176 double bondangle; // bond angle of the bond to be replaced/cut
177 double BondRescale; // rescale value for the hydrogen bond length
178 bool AllWentWell = true; // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
179 bond *FirstBond = NULL, *SecondBond = NULL; // Other bonds in double bond case to determine "other" plane
180 atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
181 int i; // loop variable
182 double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination
183 vector OrthoVector1, OrthoVector2; // temporary vectors in coordination construction
184 vector InBondVector; // vector in direction of *Bond
185 bond *Binder = NULL;
186 double *matrix;
187
188// *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
189 // create vector in direction of bond
190 InBondVector.CopyVector(&TopReplacement->x);
191 InBondVector.SubtractVector(&TopOrigin->x);
192 bondlength = InBondVector.Norm();
193
194 // is greater than typical bond distance? Then we have to correct periodically
195 // the problem is not the H being out of the box, but InBondVector have the wrong direction
196 // due to TopReplacement or Origin being on the wrong side!
197 if (bondlength > BondDistance) {
198// *out << Verbose(4) << "InBondVector is: ";
199// InBondVector.Output(out);
200// *out << endl;
201 OrthoVector1.Zero();
202 for (int i=0;i<NDIM;i++) {
203 l = TopReplacement->x.x[i] - TopOrigin->x.x[i];
204 if (fabs(l) > BondDistance) { // is component greater than bond distance
205 OrthoVector1.x[i] = (l < 0) ? -1. : +1.;
206 } // (signs are correct, was tested!)
207 }
208 matrix = ReturnFullMatrixforSymmetric(cell_size);
209 OrthoVector1.MatrixMultiplication(matrix);
210 InBondVector.SubtractVector(&OrthoVector1); // subtract just the additional translation
211 Free((void **)&matrix, "molecule::AddHydrogenReplacementAtom: *matrix");
212 bondlength = InBondVector.Norm();
213// *out << Verbose(4) << "Corrected InBondVector is now: ";
214// InBondVector.Output(out);
215// *out << endl;
216 } // periodic correction finished
217
218 InBondVector.Normalize();
219 // get typical bond length and store as scale factor for later
220 BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
221 if (BondRescale == -1) {
222 cerr << Verbose(3) << "WARNING: There is no typical bond distance for bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
223 BondRescale = bondlength;
224 } else {
225 if (!IsAngstroem)
226 BondRescale /= (1.*AtomicLengthToAngstroem);
227 }
228
229 // discern single, double and triple bonds
230 switch(TopBond->BondDegree) {
231 case 1:
232 FirstOtherAtom = new atom(); // new atom
233 FirstOtherAtom->type = elemente->FindElement(1); // element is Hydrogen
234 FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
235 FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
236 if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
237 FirstOtherAtom->father = TopReplacement;
238 BondRescale = bondlength;
239 } else {
240 FirstOtherAtom->father = NULL; // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
241 }
242 InBondVector.Scale(&BondRescale); // rescale the distance vector to Hydrogen bond length
243 FirstOtherAtom->x.CopyVector(&TopOrigin->x); // set coordination to origin ...
244 FirstOtherAtom->x.AddVector(&InBondVector); // ... and add distance vector to replacement atom
245 AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
246// *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
247// FirstOtherAtom->x.Output(out);
248// *out << endl;
249 Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
250 Binder->Cyclic = false;
251 Binder->Type = TreeEdge;
252 break;
253 case 2:
254 // determine two other bonds (warning if there are more than two other) plus valence sanity check
255 for (i=0;i<NumBond;i++) {
256 if (BondList[i] != TopBond) {
257 if (FirstBond == NULL) {
258 FirstBond = BondList[i];
259 FirstOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
260 } else if (SecondBond == NULL) {
261 SecondBond = BondList[i];
262 SecondOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
263 } else {
264 *out << Verbose(3) << "WARNING: Detected more than four bonds for atom " << TopOrigin->Name;
265 }
266 }
267 }
268 if (SecondOtherAtom == NULL) { // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
269 SecondBond = TopBond;
270 SecondOtherAtom = TopReplacement;
271 }
272 if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
273// *out << Verbose(3) << "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane." << endl;
274
275 // determine the plane of these two with the *origin
276 AllWentWell = AllWentWell && OrthoVector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
277 } else {
278 OrthoVector1.GetOneNormalVector(&InBondVector);
279 }
280 //*out << Verbose(3)<< "Orthovector1: ";
281 //OrthoVector1.Output(out);
282 //*out << endl;
283 // orthogonal vector and bond vector between origin and replacement form the new plane
284 OrthoVector1.MakeNormalVector(&InBondVector);
285 OrthoVector1.Normalize();
286 //*out << Verbose(3) << "ReScaleCheck: " << OrthoVector1.Norm() << " and " << InBondVector.Norm() << "." << endl;
287
288 // create the two Hydrogens ...
289 FirstOtherAtom = new atom();
290 SecondOtherAtom = new atom();
291 FirstOtherAtom->type = elemente->FindElement(1);
292 SecondOtherAtom->type = elemente->FindElement(1);
293 FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
294 FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
295 SecondOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
296 SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
297 FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father
298 SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father
299 bondangle = TopOrigin->type->HBondAngle[1];
300 if (bondangle == -1) {
301 *out << Verbose(3) << "WARNING: There is no typical bond angle for bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
302 bondangle = 0;
303 }
304 bondangle *= M_PI/180./2.;
305// *out << Verbose(3) << "ReScaleCheck: InBondVector ";
306// InBondVector.Output(out);
307// *out << endl;
308// *out << Verbose(3) << "ReScaleCheck: Orthovector ";
309// OrthoVector1.Output(out);
310// *out << endl;
311// *out << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
312 FirstOtherAtom->x.Zero();
313 SecondOtherAtom->x.Zero();
314 for(i=0;i<NDIM;i++) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
315 FirstOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (sin(bondangle));
316 SecondOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (-sin(bondangle));
317 }
318 FirstOtherAtom->x.Scale(&BondRescale); // rescale by correct BondDistance
319 SecondOtherAtom->x.Scale(&BondRescale);
320 //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
321 for(i=0;i<NDIM;i++) { // and make relative to origin atom
322 FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
323 SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
324 }
325 // ... and add to molecule
326 AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
327 AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
328// *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
329// FirstOtherAtom->x.Output(out);
330// *out << endl;
331// *out << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
332// SecondOtherAtom->x.Output(out);
333// *out << endl;
334 Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
335 Binder->Cyclic = false;
336 Binder->Type = TreeEdge;
337 Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
338 Binder->Cyclic = false;
339 Binder->Type = TreeEdge;
340 break;
341 case 3:
342 // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
343 FirstOtherAtom = new atom();
344 SecondOtherAtom = new atom();
345 ThirdOtherAtom = new atom();
346 FirstOtherAtom->type = elemente->FindElement(1);
347 SecondOtherAtom->type = elemente->FindElement(1);
348 ThirdOtherAtom->type = elemente->FindElement(1);
349 FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
350 FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
351 SecondOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
352 SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
353 ThirdOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
354 ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
355 FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father
356 SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father
357 ThirdOtherAtom->father = NULL; // we are just an added hydrogen with no father
358
359 // we need to vectors orthonormal the InBondVector
360 AllWentWell = AllWentWell && OrthoVector1.GetOneNormalVector(&InBondVector);
361// *out << Verbose(3) << "Orthovector1: ";
362// OrthoVector1.Output(out);
363// *out << endl;
364 AllWentWell = AllWentWell && OrthoVector2.MakeNormalVector(&InBondVector, &OrthoVector1);
365// *out << Verbose(3) << "Orthovector2: ";
366// OrthoVector2.Output(out);
367// *out << endl;
368
369 // create correct coordination for the three atoms
370 alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.; // retrieve triple bond angle from database
371 l = BondRescale; // desired bond length
372 b = 2.*l*sin(alpha); // base length of isosceles triangle
373 d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector
374 f = b/sqrt(3.); // length for OrthVector1
375 g = b/2.; // length for OrthVector2
376// *out << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
377// *out << Verbose(3) << "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << endl;
378 factors[0] = d;
379 factors[1] = f;
380 factors[2] = 0.;
381 FirstOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
382 factors[1] = -0.5*f;
383 factors[2] = g;
384 SecondOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
385 factors[2] = -g;
386 ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
387
388 // rescale each to correct BondDistance
389// FirstOtherAtom->x.Scale(&BondRescale);
390// SecondOtherAtom->x.Scale(&BondRescale);
391// ThirdOtherAtom->x.Scale(&BondRescale);
392
393 // and relative to *origin atom
394 FirstOtherAtom->x.AddVector(&TopOrigin->x);
395 SecondOtherAtom->x.AddVector(&TopOrigin->x);
396 ThirdOtherAtom->x.AddVector(&TopOrigin->x);
397
398 // ... and add to molecule
399 AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
400 AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
401 AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
402// *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
403// FirstOtherAtom->x.Output(out);
404// *out << endl;
405// *out << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
406// SecondOtherAtom->x.Output(out);
407// *out << endl;
408// *out << Verbose(4) << "Added " << *ThirdOtherAtom << " at: ";
409// ThirdOtherAtom->x.Output(out);
410// *out << endl;
411 Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
412 Binder->Cyclic = false;
413 Binder->Type = TreeEdge;
414 Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
415 Binder->Cyclic = false;
416 Binder->Type = TreeEdge;
417 Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
418 Binder->Cyclic = false;
419 Binder->Type = TreeEdge;
420 break;
421 default:
422 cerr << "ERROR: BondDegree does not state single, double or triple bond!" << endl;
423 AllWentWell = false;
424 break;
425 }
426
427// *out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
428 return AllWentWell;
429};
430
431/** Adds given atom \a *pointer from molecule list.
432 * Increases molecule::last_atom and gives last number to added atom.
433 * \param filename name and path of xyz file
434 * \return true - succeeded, false - file not found
435 */
436bool molecule::AddXYZFile(string filename)
437{
438 istringstream *input = NULL;
439 int NumberOfAtoms = 0; // atom number in xyz read
440 int i, j; // loop variables
441 atom *first = NULL; // pointer to added atom
442 char shorthand[3]; // shorthand for atom name
443 ifstream xyzfile; // xyz file
444 string line; // currently parsed line
445 double x[3]; // atom coordinates
446
447 xyzfile.open(filename.c_str());
448 if (!xyzfile)
449 return false;
450
451 getline(xyzfile,line,'\n'); // Read numer of atoms in file
452 input = new istringstream(line);
453 *input >> NumberOfAtoms;
454 cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
455 getline(xyzfile,line,'\n'); // Read comment
456 cout << Verbose(1) << "Comment: " << line << endl;
457
458 for(i=0;i<NumberOfAtoms;i++){
459 first = new atom;
460 getline(xyzfile,line,'\n');
461 istringstream *item = new istringstream(line);
462 //istringstream input(line);
463 cout << Verbose(1) << "Reading: " << line << endl;
464 *item >> shorthand;
465 *item >> x[0];
466 *item >> x[1];
467 *item >> x[2];
468 first->type = elemente->FindElement(shorthand);
469 if (first->type == NULL) {
470 cerr << "Could not parse the element at line: '" << line << "', setting to H.";
471 first->type = elemente->FindElement(1);
472 }
473 for(j=0;j<NDIM;j++)
474 first->x.x[j] = x[j];
475 AddAtom(first); // add to molecule
476 delete(item);
477 }
478 xyzfile.close();
479 delete(input);
480 return true;
481};
482
483/** Creates a copy of this molecule.
484 * \return copy of molecule
485 */
486molecule *molecule::CopyMolecule()
487{
488 molecule *copy = new molecule(elemente);
489 atom *CurrentAtom = NULL;
490 atom *LeftAtom = NULL, *RightAtom = NULL;
491 atom *Walker = NULL;
492
493 // copy all atoms
494 Walker = start;
495 while(Walker->next != end) {
496 Walker = Walker->next;
497 CurrentAtom = copy->AddCopyAtom(Walker);
498 }
499
500 // copy all bonds
501 bond *Binder = first;
502 bond *NewBond = NULL;
503 while(Binder->next != last) {
504 Binder = Binder->next;
505 // get the pendant atoms of current bond in the copy molecule
506 LeftAtom = copy->start;
507 while (LeftAtom->next != copy->end) {
508 LeftAtom = LeftAtom->next;
509 if (LeftAtom->father == Binder->leftatom)
510 break;
511 }
512 RightAtom = copy->start;
513 while (RightAtom->next != copy->end) {
514 RightAtom = RightAtom->next;
515 if (RightAtom->father == Binder->rightatom)
516 break;
517 }
518 NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
519 NewBond->Cyclic = Binder->Cyclic;
520 if (Binder->Cyclic)
521 copy->NoCyclicBonds++;
522 NewBond->Type = Binder->Type;
523 }
524 // correct fathers
525 Walker = copy->start;
526 while(Walker->next != copy->end) {
527 Walker = Walker->next;
528 if (Walker->father->father == Walker->father) // same atom in copy's father points to itself
529 Walker->father = Walker; // set father to itself (copy of a whole molecule)
530 else
531 Walker->father = Walker->father->father; // set father to original's father
532 }
533 // copy values
534 copy->CountAtoms((ofstream *)&cout);
535 copy->CountElements();
536 if (first->next != last) { // if adjaceny list is present
537 copy->BondDistance = BondDistance;
538 copy->CreateListOfBondsPerAtom((ofstream *)&cout);
539 }
540
541 return copy;
542};
543
544/** Adds a bond to a the molecule specified by two atoms, \a *first and \a *second.
545 * Also updates molecule::BondCount and molecule::NoNonBonds.
546 * \param *first first atom in bond
547 * \param *second atom in bond
548 * \return pointer to bond or NULL on failure
549 */
550bond * molecule::AddBond(atom *atom1, atom *atom2, int degree=1)
551{
552 bond *Binder = NULL;
553 if ((atom1 != NULL) && (FindAtom(atom1->nr) != NULL) && (atom2 != NULL) && (FindAtom(atom2->nr) != NULL)) {
554 Binder = new bond(atom1, atom2, degree, BondCount++);
555 if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
556 NoNonBonds++;
557 add(Binder, last);
558 } else {
559 cerr << Verbose(1) << "ERROR: Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
560 }
561 return Binder;
562};
563
564/** Remove bond from bond chain list.
565 * \todo Function not implemented yet
566 * \param *pointer bond pointer
567 * \return true - bound found and removed, false - bond not found/removed
568 */
569bool molecule::RemoveBond(bond *pointer)
570{
571 //cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
572 removewithoutcheck(pointer);
573 return true;
574};
575
576/** Remove every bond from bond chain list that atom \a *BondPartner is a constituent of.
577 * \todo Function not implemented yet
578 * \param *BondPartner atom to be removed
579 * \return true - bounds found and removed, false - bonds not found/removed
580 */
581bool molecule::RemoveBonds(atom *BondPartner)
582{
583 cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
584 return false;
585};
586
587/** Sets the molecule::cell_size to the components of \a *dim (rectangular box)
588 * \param *dim vector class
589 */
590void molecule::SetBoxDimension(vector *dim)
591{
592 cell_size[0] = dim->x[0];
593 cell_size[1] = 0.;
594 cell_size[2] = dim->x[1];
595 cell_size[3] = 0.;
596 cell_size[4] = 0.;
597 cell_size[5] = dim->x[2];
598};
599
600/** Centers the edge of the atoms at (0,0,0).
601 * \param *out output stream for debugging
602 * \param *max coordinates of other edge, specifying box dimensions.
603 */
604void molecule::CenterEdge(ofstream *out, vector *max)
605{
606 vector *min = new vector;
607
608// *out << Verbose(3) << "Begin of CenterEdge." << endl;
609 atom *ptr = start->next; // start at first in list
610 if (ptr != end) { //list not empty?
611 for (int i=0;i<NDIM;i++) {
612 max->x[i] = ptr->x.x[i];
613 min->x[i] = ptr->x.x[i];
614 }
615 while (ptr->next != end) { // continue with second if present
616 ptr = ptr->next;
617 //ptr->Output(1,1,out);
618 for (int i=0;i<NDIM;i++) {
619 max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
620 min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i];
621 }
622 }
623// *out << Verbose(4) << "Maximum is ";
624// max->Output(out);
625// *out << ", Minimum is ";
626// min->Output(out);
627// *out << endl;
628
629 for (int i=0;i<NDIM;i++) {
630 min->x[i] *= -1.;
631 max->x[i] += min->x[i];
632 }
633 Translate(min);
634 }
635 delete(min);
636// *out << Verbose(3) << "End of CenterEdge." << endl;
637};
638
639/** Centers the center of the atoms at (0,0,0).
640 * \param *out output stream for debugging
641 * \param *center return vector for translation vector
642 */
643void molecule::CenterOrigin(ofstream *out, vector *center)
644{
645 int Num = 0;
646 atom *ptr = start->next; // start at first in list
647
648 for(int i=0;i<NDIM;i++) // zero center vector
649 center->x[i] = 0.;
650
651 if (ptr != end) { //list not empty?
652 while (ptr->next != end) { // continue with second if present
653 ptr = ptr->next;
654 Num++;
655 center->AddVector(&ptr->x);
656 }
657 center->Scale(-1./Num); // divide through total number (and sign for direction)
658 Translate(center);
659 }
660};
661
662/** Centers the center of gravity of the atoms at (0,0,0).
663 * \param *out output stream for debugging
664 * \param *center return vector for translation vector
665 */
666void molecule::CenterGravity(ofstream *out, vector *center)
667{
668 double Num = 0;
669 atom *ptr = start->next; // start at first in list
670 vector tmp;
671
672 for(int i=0;i<NDIM;i++) // zero center vector
673 center->x[i] = 0.;
674
675 if (ptr != end) { //list not empty?
676 while (ptr->next != end) { // continue with second if present
677 ptr = ptr->next;
678 Num += ptr->type->mass;
679 tmp.CopyVector(&ptr->x);
680 tmp.Scale(ptr->type->mass); // scale by mass
681 center->AddVector(&tmp);
682 }
683 center->Scale(-1./Num); // divide through total mass (and sign for direction)
684 Translate(center);
685 }
686};
687
688/** Scales all atoms by \a *factor.
689 * \param *factor pointer to scaling factor
690 */
691void molecule::Scale(double **factor)
692{
693 atom *ptr = start;
694
695 while (ptr->next != end) {
696 ptr = ptr->next;
697 ptr->x.Scale(factor);
698 }
699};
700
701/** Translate all atoms by given vector.
702 * \param trans[] translation vector.
703 */
704void molecule::Translate(const vector *trans)
705{
706 atom *ptr = start;
707
708 while (ptr->next != end) {
709 ptr = ptr->next;
710 ptr->x.Translate(trans);
711 }
712};
713
714/** Mirrors all atoms against a given plane.
715 * \param n[] normal vector of mirror plane.
716 */
717void molecule::Mirror(const vector *n)
718{
719 atom *ptr = start;
720
721 while (ptr->next != end) {
722 ptr = ptr->next;
723 ptr->x.Mirror(n);
724 }
725};
726
727/** Determines center of gravity (yet not considering atom masses).
728 * \param CenterOfGravity reference to return vector
729 */
730void molecule::DetermineCenterOfGravity(vector &CenterOfGravity)
731{
732 atom *Walker = start;
733 bond *Binder = NULL;
734 double *matrix = ReturnFullMatrixforSymmetric(cell_size);
735 double tmp;
736 bool flag;
737 vector TestVector, TranslationVector;
738
739 do {
740 CenterOfGravity.Zero();
741 flag = true;
742 while (Walker->next != end) {
743 Walker = Walker->next;
744#ifdef ADDHYDROGEN
745 if (Walker->type->Z != 1) {
746#endif
747 TestVector.CopyVector(&Walker->x);
748 TestVector.InverseMatrixMultiplication(matrix);
749 TranslationVector.Zero();
750 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
751 Binder = ListOfBondsPerAtom[Walker->nr][i];
752 if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
753 for (int j=0;j<NDIM;j++) {
754 tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
755 if ((fabs(tmp)) > BondDistance) {
756 flag = false;
757 cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl;
758 if (tmp > 0)
759 TranslationVector.x[j] -= 1.;
760 else
761 TranslationVector.x[j] += 1.;
762 }
763 }
764 }
765 TestVector.AddVector(&TranslationVector);
766 TestVector.MatrixMultiplication(matrix);
767 CenterOfGravity.AddVector(&TestVector);
768 cout << Verbose(1) << "Vector is: ";
769 TestVector.Output((ofstream *)&cout);
770 cout << endl;
771#ifdef ADDHYDROGEN
772 // now also change all hydrogens
773 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
774 Binder = ListOfBondsPerAtom[Walker->nr][i];
775 if (Binder->GetOtherAtom(Walker)->type->Z == 1) {
776 TestVector.CopyVector(&Binder->GetOtherAtom(Walker)->x);
777 TestVector.InverseMatrixMultiplication(matrix);
778 TestVector.AddVector(&TranslationVector);
779 TestVector.MatrixMultiplication(matrix);
780 CenterOfGravity.AddVector(&TestVector);
781 cout << Verbose(1) << "Hydrogen Vector is: ";
782 TestVector.Output((ofstream *)&cout);
783 cout << endl;
784 }
785 }
786 }
787#endif
788 }
789 } while (!flag);
790 Free((void **)&matrix, "molecule::DetermineCenterOfGravity: *matrix");
791 CenterOfGravity.Scale(1./(double)AtomCount);
792};
793
794/** Align all atoms in such a manner that given vector \a *n is along z axis.
795 * \param n[] alignment vector.
796 */
797void molecule::Align(vector *n)
798{
799 atom *ptr = start;
800 double alpha, tmp;
801 vector z_axis;
802 z_axis.x[0] = 0.;
803 z_axis.x[1] = 0.;
804 z_axis.x[2] = 1.;
805
806 // rotate on z-x plane
807 cout << Verbose(0) << "Begin of Aligning all atoms." << endl;
808 alpha = atan(-n->x[0]/n->x[2]);
809 cout << Verbose(1) << "Z-X-angle: " << alpha << " ... ";
810 while (ptr->next != end) {
811 ptr = ptr->next;
812 tmp = ptr->x.x[0];
813 ptr->x.x[0] = cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
814 ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
815 }
816 // rotate n vector
817 tmp = n->x[0];
818 n->x[0] = cos(alpha) * tmp + sin(alpha) * n->x[2];
819 n->x[2] = -sin(alpha) * tmp + cos(alpha) * n->x[2];
820 cout << Verbose(1) << "alignment vector after first rotation: ";
821 n->Output((ofstream *)&cout);
822 cout << endl;
823
824 // rotate on z-y plane
825 ptr = start;
826 alpha = atan(-n->x[1]/n->x[2]);
827 cout << Verbose(1) << "Z-Y-angle: " << alpha << " ... ";
828 while (ptr->next != end) {
829 ptr = ptr->next;
830 tmp = ptr->x.x[1];
831 ptr->x.x[1] = cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
832 ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
833 }
834 // rotate n vector (for consistency check)
835 tmp = n->x[1];
836 n->x[1] = cos(alpha) * tmp + sin(alpha) * n->x[2];
837 n->x[2] = -sin(alpha) * tmp + cos(alpha) * n->x[2];
838
839 cout << Verbose(1) << "alignment vector after second rotation: ";
840 n->Output((ofstream *)&cout);
841 cout << Verbose(1) << endl;
842 cout << Verbose(0) << "End of Aligning all atoms." << endl;
843};
844
845/** Removes atom from molecule list.
846 * \param *pointer atom to be removed
847 * \return true - succeeded, false - atom not found in list
848 */
849bool molecule::RemoveAtom(atom *pointer)
850{
851 if (ElementsInMolecule[pointer->type->Z] != 0) // this would indicate an error
852 ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
853 else
854 cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
855 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element?
856 ElementCount--;
857 return remove(pointer, start, end);
858};
859
860/** Removes every atom from molecule list.
861 * \return true - succeeded, false - atom not found in list
862 */
863bool molecule::CleanupMolecule()
864{
865 return (cleanup(start,end) && cleanup(first,last));
866};
867
868/** Finds an atom specified by its continuous number.
869 * \param Nr number of atom withim molecule
870 * \return pointer to atom or NULL
871 */
872atom * molecule::FindAtom(int Nr) const{
873 atom * walker = find(&Nr, start,end);
874 if (walker != NULL) {
875 //cout << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
876 return walker;
877 } else {
878 cout << Verbose(0) << "Atom not found in list." << endl;
879 return NULL;
880 }
881};
882
883/** Asks for atom number, and checks whether in list.
884 * \param *text question before entering
885 */
886atom * molecule::AskAtom(char *text)
887{
888 int No;
889 atom *ion = NULL;
890 do {
891 //cout << Verbose(0) << "============Atom list==========================" << endl;
892 //mol->Output((ofstream *)&cout);
893 //cout << Verbose(0) << "===============================================" << endl;
894 cout << Verbose(0) << text;
895 cin >> No;
896 ion = this->FindAtom(No);
897 } while (ion == NULL);
898 return ion;
899};
900
901/** Checks if given coordinates are within cell volume.
902 * \param *x array of coordinates
903 * \return true - is within, false - out of cell
904 */
905bool molecule::CheckBounds(const vector *x) const
906{
907 bool result = true;
908 int j =-1;
909 for (int i=0;i<3;i++) {
910 j += i+1;
911 result = result && ((x->x[i] >= 0) && (x->x[i] < cell_size[j]));
912 }
913 //return result;
914 return true; /// probably not gonna use the check no more
915};
916
917/** Calculates sum over least square distance to line hidden in \a *x.
918 * \param *x offset and direction vector
919 * \param *params pointer to lsq_params structure
920 * \return \f$ sum_i^N | y_i - (a + t_i b)|^2\f$
921 */
922double LeastSquareDistance (const gsl_vector * x, void * params)
923{
924 double res = 0, t;
925 vector a,b,c,d;
926 struct lsq_params *par = (struct lsq_params *)params;
927 atom *ptr = par->mol->start;
928
929 // initialize vectors
930 a.x[0] = gsl_vector_get(x,0);
931 a.x[1] = gsl_vector_get(x,1);
932 a.x[2] = gsl_vector_get(x,2);
933 b.x[0] = gsl_vector_get(x,3);
934 b.x[1] = gsl_vector_get(x,4);
935 b.x[2] = gsl_vector_get(x,5);
936 // go through all atoms
937 while (ptr != par->mol->end) {
938 ptr = ptr->next;
939 if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
940 c.CopyVector(&ptr->x); // copy vector to temporary one
941 c.SubtractVector(&a); // subtract offset vector
942 t = c.ScalarProduct(&b); // get direction parameter
943 d.CopyVector(&b); // and create vector
944 d.Scale(&t);
945 c.SubtractVector(&d); // ... yielding distance vector
946 res += d.ScalarProduct((const vector *)&d); // add squared distance
947 }
948 }
949 return res;
950};
951
952/** By minimizing the least square distance gains alignment vector.
953 * \bug this is not yet working properly it seems
954 */
955void molecule::GetAlignVector(struct lsq_params * par) const
956{
957 int np = 6;
958
959 const gsl_multimin_fminimizer_type *T =
960 gsl_multimin_fminimizer_nmsimplex;
961 gsl_multimin_fminimizer *s = NULL;
962 gsl_vector *ss;
963 gsl_multimin_function minex_func;
964
965 size_t iter = 0, i;
966 int status;
967 double size;
968
969 /* Initial vertex size vector */
970 ss = gsl_vector_alloc (np);
971
972 /* Set all step sizes to 1 */
973 gsl_vector_set_all (ss, 1.0);
974
975 /* Starting point */
976 par->x = gsl_vector_alloc (np);
977 par->mol = this;
978
979 gsl_vector_set (par->x, 0, 0.0); // offset
980 gsl_vector_set (par->x, 1, 0.0);
981 gsl_vector_set (par->x, 2, 0.0);
982 gsl_vector_set (par->x, 3, 0.0); // direction
983 gsl_vector_set (par->x, 4, 0.0);
984 gsl_vector_set (par->x, 5, 1.0);
985
986 /* Initialize method and iterate */
987 minex_func.f = &LeastSquareDistance;
988 minex_func.n = np;
989 minex_func.params = (void *)par;
990
991 s = gsl_multimin_fminimizer_alloc (T, np);
992 gsl_multimin_fminimizer_set (s, &minex_func, par->x, ss);
993
994 do
995 {
996 iter++;
997 status = gsl_multimin_fminimizer_iterate(s);
998
999 if (status)
1000 break;
1001
1002 size = gsl_multimin_fminimizer_size (s);
1003 status = gsl_multimin_test_size (size, 1e-2);
1004
1005 if (status == GSL_SUCCESS)
1006 {
1007 printf ("converged to minimum at\n");
1008 }
1009
1010 printf ("%5d ", (int)iter);
1011 for (i = 0; i < (size_t)np; i++)
1012 {
1013 printf ("%10.3e ", gsl_vector_get (s->x, i));
1014 }
1015 printf ("f() = %7.3f size = %.3f\n", s->fval, size);
1016 }
1017 while (status == GSL_CONTINUE && iter < 100);
1018
1019 for (i=0;i<(size_t)np;i++)
1020 gsl_vector_set(par->x, i, gsl_vector_get(s->x, i));
1021 //gsl_vector_free(par->x);
1022 gsl_vector_free(ss);
1023 gsl_multimin_fminimizer_free (s);
1024};
1025
1026/** Prints molecule to *out.
1027 * \param *out output stream
1028 */
1029bool molecule::Output(ofstream *out)
1030{
1031 element *runner = elemente->start;
1032 atom *walker = NULL;
1033 int ElementNo, AtomNo;
1034 CountElements();
1035
1036 if (out == NULL) {
1037 return false;
1038 } else {
1039 *out << "#Ion_TypeNr._Nr.R[0] R[1] R[2] MoveType (0 MoveIon, 1 FixedIon)" << endl;
1040 ElementNo = 0;
1041 while (runner->next != elemente->end) { // go through every element
1042 runner = runner->next;
1043 if (ElementsInMolecule[runner->Z]) { // if this element got atoms
1044 ElementNo++;
1045 AtomNo = 0;
1046 walker = start;
1047 while (walker->next != end) { // go through every atom of this element
1048 walker = walker->next;
1049 if (walker->type == runner) { // if this atom fits to element
1050 AtomNo++;
1051 walker->Output(ElementNo, AtomNo, out);
1052 }
1053 }
1054 }
1055 }
1056 return true;
1057 }
1058};
1059
1060/** Outputs contents of molecule::ListOfBondsPerAtom.
1061 * \param *out output stream
1062 */
1063void molecule::OutputListOfBonds(ofstream *out) const
1064{
1065 *out << Verbose(2) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
1066 atom *Walker = start;
1067 while (Walker->next != end) {
1068 Walker = Walker->next;
1069#ifdef ADDHYDROGEN
1070 if (Walker->type->Z != 1) { // regard only non-hydrogen
1071#endif
1072 *out << Verbose(2) << "Atom " << Walker->Name << " has Bonds: "<<endl;
1073 for(int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {
1074 *out << Verbose(3) << *(ListOfBondsPerAtom)[Walker->nr][j] << endl;
1075 }
1076#ifdef ADDHYDROGEN
1077 }
1078#endif
1079 }
1080 *out << endl;
1081};
1082
1083/** Output of element before the actual coordination list.
1084 * \param *out stream pointer
1085 */
1086bool molecule::Checkout(ofstream *out) const
1087{
1088 return elemente->Checkout(out, ElementsInMolecule);
1089};
1090
1091/** Prints molecule to *out as xyz file.
1092 * \param *out output stream
1093 */
1094bool molecule::OutputXYZ(ofstream *out) const
1095{
1096 atom *walker = NULL;
1097 int No = 0;
1098 time_t now;
1099
1100 now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time'
1101 walker = start;
1102 while (walker->next != end) { // go through every atom and count
1103 walker = walker->next;
1104 No++;
1105 }
1106 if (out != NULL) {
1107 *out << No << "\n\tCreated by molecuilder on " << ctime(&now);
1108 walker = start;
1109 while (walker->next != end) { // go through every atom of this element
1110 walker = walker->next;
1111 walker->OutputXYZLine(out);
1112 }
1113 return true;
1114 } else
1115 return false;
1116};
1117
1118/** Brings molecule::AtomCount and atom::*Name up-to-date.
1119 * \param *out output stream for debugging
1120 */
1121void molecule::CountAtoms(ofstream *out)
1122{
1123 int i = 0;
1124 atom *Walker = start;
1125 while (Walker->next != end) {
1126 Walker = Walker->next;
1127 i++;
1128 }
1129 if ((AtomCount == 0) || (i != AtomCount)) {
1130 *out << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
1131 AtomCount = i;
1132
1133 // count NonHydrogen atoms and give each atom a unique name
1134 if (AtomCount != 0) {
1135 i=0;
1136 NoNonHydrogen = 0;
1137 Walker = start;
1138 while (Walker->next != end) {
1139 Walker = Walker->next;
1140 Walker->nr = i; // update number in molecule (for easier referencing in FragmentMolecule lateron)
1141 if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
1142 NoNonHydrogen++;
1143 Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
1144 Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
1145 sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
1146 *out << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
1147 i++;
1148 }
1149 } else
1150 *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
1151 }
1152};
1153
1154/** Brings molecule::ElementCount and molecule::ElementsInMolecule up-to-date.
1155 */
1156void molecule::CountElements()
1157{
1158 int i = 0;
1159 for(i=0;i<MAX_ELEMENTS;i++)
1160 ElementsInMolecule[i] = 0;
1161 ElementCount = 0;
1162
1163 atom *walker = start;
1164 while (walker->next != end) {
1165 walker = walker->next;
1166 ElementsInMolecule[walker->type->Z]++;
1167 i++;
1168 }
1169 for(i=0;i<MAX_ELEMENTS;i++)
1170 ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
1171};
1172
1173/** Counts all cyclic bonds and returns their number.
1174 * \note Hydrogen bonds can never by cyclic, thus no check for that
1175 * \param *out output stream for debugging
1176 * \return number opf cyclic bonds
1177 */
1178int molecule::CountCyclicBonds(ofstream *out)
1179{
1180 int No = 0;
1181 int *MinimumRingSize = NULL;
1182 MoleculeLeafClass *Subgraphs = NULL;
1183 bond *Binder = first;
1184 if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
1185 *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
1186 Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
1187 while (Subgraphs->next != NULL) {
1188 Subgraphs = Subgraphs->next;
1189 delete(Subgraphs->previous);
1190 }
1191 delete(Subgraphs);
1192 delete[](MinimumRingSize);
1193 }
1194 while(Binder->next != last) {
1195 Binder = Binder->next;
1196 if (Binder->Cyclic)
1197 No++;
1198 }
1199 return No;
1200};
1201
1202/** Returns Shading as a char string.
1203 * \param color the Shading
1204 * \return string of the flag
1205 */
1206char * molecule::GetColor(enum Shading color)
1207{
1208 switch(color) {
1209 case white:
1210 return "white";
1211 break;
1212 case lightgray:
1213 return "lightgray";
1214 break;
1215 case darkgray:
1216 return "darkgray";
1217 break;
1218 case black:
1219 return "black";
1220 break;
1221 default:
1222 return "uncolored";
1223 break;
1224 };
1225};
1226
1227
1228/** Counts necessary number of valence electrons and returns number and SpinType.
1229 * \param configuration containing everything
1230 */
1231void molecule::CalculateOrbitals(class config &configuration)
1232{
1233 configuration.MaxPsiDouble = configuration.PsiMaxNoDown = configuration.PsiMaxNoUp = configuration.PsiType = 0;
1234 for(int i=0;i<MAX_ELEMENTS;i++) {
1235 if (ElementsInMolecule[i] != 0) {
1236 //cout << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
1237 configuration.MaxPsiDouble += ElementsInMolecule[i]*((int)elemente->FindElement(i)->Valence);
1238 }
1239 }
1240 configuration.PsiMaxNoDown = configuration.MaxPsiDouble/2 + (configuration.MaxPsiDouble % 2);
1241 configuration.PsiMaxNoUp = configuration.MaxPsiDouble/2;
1242 configuration.MaxPsiDouble /= 2;
1243 configuration.PsiType = (configuration.PsiMaxNoDown == configuration.PsiMaxNoUp) ? 0 : 1;
1244 if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2)) {
1245 configuration.ProcPEGamma /= 2;
1246 configuration.ProcPEPsi *= 2;
1247 } else {
1248 configuration.ProcPEGamma *= configuration.ProcPEPsi;
1249 configuration.ProcPEPsi = 1;
1250 }
1251 configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
1252};
1253
1254/** Creates an adjacency list of the molecule.
1255 * Generally, we use the CSD approach to bond recognition, that is the the distance
1256 * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
1257 * a threshold t = 0.4 Angstroem.
1258 * To make it O(N log N) the function uses the linked-cell technique as follows:
1259 * The procedure is step-wise:
1260 * -# Remove every bond in list
1261 * -# Count the atoms in the molecule with CountAtoms()
1262 * -# partition cell into smaller linked cells of size \a bonddistance
1263 * -# put each atom into its corresponding cell
1264 * -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
1265 * -# create the list of bonds via CreateListOfBondsPerAtom()
1266 * -# correct the bond degree iteratively (single->double->triple bond)
1267 * -# finally print the bond list to \a *out if desired
1268 * \param *out out stream for printing the matrix, NULL if no output
1269 * \param bonddistance length of linked cells (i.e. maximum minimal length checked)
1270 */
1271void molecule::CreateAdjacencyList(ofstream *out, double bonddistance)
1272{
1273 atom *Walker = NULL, *OtherWalker = NULL;
1274 int No, NoBonds;
1275 int NumberCells, divisor[NDIM], n[NDIM], N[NDIM], index, Index, j;
1276 molecule **CellList;
1277 double distance, MinDistance, MaxDistance;
1278 double *matrix = ReturnFullMatrixforSymmetric(cell_size);
1279 vector x;
1280
1281 BondDistance = bonddistance;
1282 *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
1283 // remove every bond from the list
1284 if ((first->next != last) && (last->previous != first)) { // there are bonds present
1285 cleanup(first,last);
1286 }
1287
1288 // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
1289 CountAtoms(out);
1290 *out << Verbose(1) << "AtomCount " << AtomCount << "." << endl;
1291
1292 if (AtomCount != 0) {
1293 // 1. find divisor for each axis, such that a sphere with radius of at least bonddistance can be placed into each cell
1294 j=-1;
1295 for (int i=0;i<NDIM;i++) {
1296 j += i+1;
1297 divisor[i] = (int)floor(cell_size[j]/bonddistance); // take smaller value such that size of linked cell is at least bonddistance
1298 *out << Verbose(1) << "divisor[" << i << "] = " << divisor[i] << "." << endl;
1299 }
1300 // 2a. allocate memory for the cell list
1301 NumberCells = divisor[0]*divisor[1]*divisor[2];
1302 *out << Verbose(1) << "Allocating " << NumberCells << " cells." << endl;
1303 CellList = (molecule **) Malloc(sizeof(molecule *)*NumberCells, "molecule::CreateAdjacencyList - ** CellList");
1304 for (int i=0;i<NumberCells;i++)
1305 CellList[i] = NULL;
1306
1307 // 2b. put all atoms into its corresponding list
1308 Walker = start;
1309 while(Walker->next != end) {
1310 Walker = Walker->next;
1311 //*out << Verbose(1) << "Current atom is " << *Walker << " with coordinates ";
1312 //Walker->x.Output(out);
1313 //*out << "." << endl;
1314 // compute the cell by the atom's coordinates
1315 j=-1;
1316 for (int i=0;i<NDIM;i++) {
1317 j += i+1;
1318 x.CopyVector(&(Walker->x));
1319 x.KeepPeriodic(out, matrix);
1320 n[i] = (int)floor(x.x[i]/cell_size[j]*(double)divisor[i]);
1321 }
1322 index = n[2] + (n[1] + n[0] * divisor[1]) * divisor[2];
1323 *out << Verbose(1) << "Atom " << *Walker << " goes into cell number [" << n[0] << "," << n[1] << "," << n[2] << "] = " << index << "." << endl;
1324 // add copy atom to this cell
1325 if (CellList[index] == NULL) // allocate molecule if not done
1326 CellList[index] = new molecule(elemente);
1327 OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
1328 //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
1329 }
1330 //for (int i=0;i<NumberCells;i++)
1331 //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
1332
1333 // 3a. go through every cell
1334 for (N[0]=0;N[0]<divisor[0];N[0]++)
1335 for (N[1]=0;N[1]<divisor[1];N[1]++)
1336 for (N[2]=0;N[2]<divisor[2];N[2]++) {
1337 Index = N[2] + (N[1] + N[0] * divisor[1]) * divisor[2];
1338 if (CellList[Index] != NULL) { // if there atoms in this cell
1339 //*out << Verbose(1) << "Current cell is " << Index << "." << endl;
1340 // 3b. for every atom therein
1341 Walker = CellList[Index]->start;
1342 while (Walker->next != CellList[Index]->end) { // go through every atom
1343 Walker = Walker->next;
1344 //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
1345 // 3c. check for possible bond between each atom in this and every one in the 27 cells
1346 for (n[0]=-1;n[0]<=1;n[0]++)
1347 for (n[1]=-1;n[1]<=1;n[1]++)
1348 for (n[2]=-1;n[2]<=1;n[2]++) {
1349 // compute the index of this comparison cell and make it periodic
1350 index = ((N[2]+n[2]+divisor[2])%divisor[2]) + (((N[1]+n[1]+divisor[1])%divisor[1]) + ((N[0]+n[0]+divisor[0])%divisor[0]) * divisor[1]) * divisor[2];
1351 //*out << Verbose(1) << "Number of comparison cell is " << index << "." << endl;
1352 if (CellList[index] != NULL) { // if there are any atoms in this cell
1353 OtherWalker = CellList[index]->start;
1354 while(OtherWalker->next != CellList[index]->end) { // go through every atom in this cell
1355 OtherWalker = OtherWalker->next;
1356 //*out << Verbose(0) << "Current comparison atom is " << *OtherWalker << "." << endl;
1357 /// \todo periodic check is missing here!
1358 //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
1359 MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
1360 MaxDistance = MinDistance + BONDTHRESHOLD;
1361 MinDistance -= BONDTHRESHOLD;
1362 distance = OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size);
1363 if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
1364 *out << Verbose(0) << "Adding Bond between " << *Walker << " and " << *OtherWalker << "." << endl;
1365 AddBond(Walker->father, OtherWalker->father, 1); // also increases molecule::BondCount
1366 } else {
1367 //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
1368 }
1369 }
1370 }
1371 }
1372 }
1373 }
1374 }
1375 // 4. free the cell again
1376 for (int i=0;i<NumberCells;i++)
1377 if (CellList[i] != NULL) {
1378 delete(CellList[i]);
1379 }
1380 Free((void **)&CellList, "molecule::CreateAdjacencyList - ** CellList");
1381
1382 // create the adjacency list per atom
1383 CreateListOfBondsPerAtom(out);
1384
1385 // correct Bond degree of each bond by checking of updated(!) sum of bond degree for an atom match its valence count
1386 // bond degrres are correctled iteratively by one, so that 2-2 instead of 1-3 or 3-1 corrections are favoured: We want
1387 // a rather symmetric distribution of higher bond degrees
1388 if (BondCount != 0) {
1389 NoCyclicBonds = 0;
1390 *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
1391 do {
1392 No = 0; // No acts as breakup flag (if 1 we still continue)
1393 Walker = start;
1394 while (Walker->next != end) { // go through every atom
1395 Walker = Walker->next;
1396 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
1397 OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
1398 // count valence of first partner (updated!), might have changed during last bond partner
1399 NoBonds = 0;
1400 for(j=0;j<NumberOfBondsPerAtom[Walker->nr];j++)
1401 NoBonds += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
1402 //*out << Verbose(3) << "Walker: " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
1403 if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check NoBonds of other atom
1404 // count valence of second partner
1405 NoBonds = 0;
1406 for(j=0;j<NumberOfBondsPerAtom[OtherWalker->nr];j++)
1407 NoBonds += ListOfBondsPerAtom[OtherWalker->nr][j]->BondDegree;
1408 //*out << Verbose(3) << "OtherWalker: " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
1409 if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) // increase bond degree by one
1410 ListOfBondsPerAtom[Walker->nr][i]->BondDegree++;
1411 }
1412 }
1413 }
1414 } while (No);
1415 *out << " done." << endl;
1416 } else
1417 *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
1418 *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << "." << endl;
1419
1420// // output bonds for debugging (if bond chain list was correctly installed)
1421// *out << Verbose(1) << endl << "From contents of bond chain list:";
1422// bond *Binder = first;
1423// while(Binder->next != last) {
1424// Binder = Binder->next;
1425// *out << *Binder << "\t" << endl;
1426// }
1427// *out << endl;
1428 } else
1429 *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
1430 *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
1431 Free((void **)&matrix, "molecule::CreateAdjacencyList: *matrix");
1432};
1433
1434/** Performs a Depth-First search on this molecule.
1435 * Marks bonds in molecule as cyclic, bridge, ... and atoms as
1436 * articulations points, ...
1437 * We use the algorithm from [Even, Graph Algorithms, p.62].
1438 * \param *out output stream for debugging
1439 * \param ReturnStack true - return pointer to atom stack of separable components, false - return NULL
1440 * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found
1441 * \return list of each disconnected subgraph as an individual molecule class structure
1442 */
1443MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int *&MinimumRingSize)
1444{
1445 class StackClass<atom *> *AtomStack;
1446 AtomStack = new StackClass<atom *>(AtomCount);
1447 class StackClass<bond *> *BackEdgeStack = new StackClass<bond *> (BondCount);
1448 MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
1449 MoleculeLeafClass *LeafWalker = SubGraphs;
1450 int CurrentGraphNr = 0, OldGraphNr;
1451 int ComponentNumber = 0;
1452 atom *Walker = NULL, *OtherAtom = NULL, *Root = start->next;
1453 bond *Binder = NULL;
1454 bool BackStepping = false;
1455
1456 *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
1457
1458 ResetAllBondsToUnused();
1459 ResetAllAtomNumbers();
1460 InitComponentNumbers();
1461 BackEdgeStack->ClearStack();
1462 while (Root != end) { // if there any atoms at all
1463 // (1) mark all edges unused, empty stack, set atom->GraphNr = 0 for all
1464 AtomStack->ClearStack();
1465
1466 // put into new subgraph molecule and add this to list of subgraphs
1467 LeafWalker = new MoleculeLeafClass(LeafWalker);
1468 LeafWalker->Leaf = new molecule(elemente);
1469 LeafWalker->Leaf->AddCopyAtom(Root);
1470
1471 OldGraphNr = CurrentGraphNr;
1472 Walker = Root;
1473 do { // (10)
1474 do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom
1475 if (!BackStepping) { // if we don't just return from (8)
1476 Walker->GraphNr = CurrentGraphNr;
1477 Walker->LowpointNr = CurrentGraphNr;
1478 *out << Verbose(1) << "Setting Walker[" << Walker->Name << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "." << endl;
1479 AtomStack->Push(Walker);
1480 CurrentGraphNr++;
1481 }
1482 do { // (3) if Walker has no unused egdes, go to (5)
1483 BackStepping = false; // reset backstepping flag for (8)
1484 if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
1485 Binder = FindNextUnused(Walker);
1486 if (Binder == NULL)
1487 break;
1488 *out << Verbose(2) << "Current Unused Bond is " << *Binder << "." << endl;
1489 // (4) Mark Binder used, ...
1490 Binder->MarkUsed(black);
1491 OtherAtom = Binder->GetOtherAtom(Walker);
1492 *out << Verbose(2) << "(4) OtherAtom is " << OtherAtom->Name << "." << endl;
1493 if (OtherAtom->GraphNr != -1) {
1494 // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
1495 Binder->Type = BackEdge;
1496 BackEdgeStack->Push(Binder);
1497 Walker->LowpointNr = ( Walker->LowpointNr < OtherAtom->GraphNr ) ? Walker->LowpointNr : OtherAtom->GraphNr;
1498 *out << Verbose(3) << "(4a) Visited: Setting Lowpoint of Walker[" << Walker->Name << "] to " << Walker->LowpointNr << "." << endl;
1499 } else {
1500 // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
1501 Binder->Type = TreeEdge;
1502 OtherAtom->Ancestor = Walker;
1503 Walker = OtherAtom;
1504 *out << Verbose(3) << "(4b) Not Visited: OtherAtom[" << OtherAtom->Name << "]'s Ancestor is now " << OtherAtom->Ancestor->Name << ", Walker is OtherAtom " << OtherAtom->Name << "." << endl;
1505 break;
1506 }
1507 Binder = NULL;
1508 } while (1); // (3)
1509 if (Binder == NULL) {
1510 *out << Verbose(2) << "No more Unused Bonds." << endl;
1511 break;
1512 } else
1513 Binder = NULL;
1514 } while (1); // (2)
1515
1516 // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
1517 if ((Walker == Root) && (Binder == NULL))
1518 break;
1519
1520 // (5) if Ancestor of Walker is ...
1521 *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
1522 if (Walker->Ancestor->GraphNr != Root->GraphNr) {
1523 // (6) (Ancestor of Walker is not Root)
1524 if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
1525 // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
1526 Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
1527 *out << Verbose(2) << "(6) Setting Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "." << endl;
1528 } else {
1529 // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
1530 Walker->Ancestor->SeparationVertex = true;
1531 *out << Verbose(2) << "(7) Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "]'s is a separating vertex, creating component." << endl;
1532 SetNextComponentNumber(Walker->Ancestor, ComponentNumber);
1533 *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Ancestor's Compont is " << ComponentNumber << "." << endl;
1534 SetNextComponentNumber(Walker, ComponentNumber);
1535 *out << Verbose(3) << "(7) Walker[" << Walker->Name << "]'s Compont is " << ComponentNumber << "." << endl;
1536 do {
1537 OtherAtom = AtomStack->PopLast();
1538 LeafWalker->Leaf->AddCopyAtom(OtherAtom);
1539 SetNextComponentNumber(OtherAtom, ComponentNumber);
1540 *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << ComponentNumber << "." << endl;
1541 } while (OtherAtom != Walker);
1542 ComponentNumber++;
1543 }
1544 // (8) Walker becomes its Ancestor, go to (3)
1545 *out << Verbose(2) << "(8) Walker[" << Walker->Name << "] is now its Ancestor " << Walker->Ancestor->Name << ", backstepping. " << endl;
1546 Walker = Walker->Ancestor;
1547 BackStepping = true;
1548 }
1549 if (!BackStepping) { // coming from (8) want to go to (3)
1550 // (9) remove all from stack till Walker (including), these and Root form a component
1551 AtomStack->Output(out);
1552 SetNextComponentNumber(Root, ComponentNumber);
1553 *out << Verbose(3) << "(9) Root[" << Root->Name << "]'s Component is " << ComponentNumber << "." << endl;
1554 SetNextComponentNumber(Walker, ComponentNumber);
1555 *out << Verbose(3) << "(9) Walker[" << Walker->Name << "]'s Component is " << ComponentNumber << "." << endl;
1556 do {
1557 OtherAtom = AtomStack->PopLast();
1558 LeafWalker->Leaf->AddCopyAtom(OtherAtom);
1559 SetNextComponentNumber(OtherAtom, ComponentNumber);
1560 *out << Verbose(3) << "(7) Other[" << OtherAtom->Name << "]'s Compont is " << ComponentNumber << "." << endl;
1561 } while (OtherAtom != Walker);
1562 ComponentNumber++;
1563
1564 // (11) Root is separation vertex, set Walker to Root and go to (4)
1565 Walker = Root;
1566 Binder = FindNextUnused(Walker);
1567 *out << Verbose(1) << "(10) Walker is Root[" << Root->Name << "], next Unused Bond is " << Binder << "." << endl;
1568 if (Binder != NULL) { // Root is separation vertex
1569 *out << Verbose(1) << "(11) Root is a separation vertex." << endl;
1570 Walker->SeparationVertex = true;
1571 }
1572 }
1573 } while ((BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
1574
1575 // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
1576 *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
1577 LeafWalker->Leaf->Output(out);
1578 *out << endl;
1579
1580 // step on to next root
1581 while ((Root != end) && (Root->GraphNr != -1)) {
1582 //*out << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
1583 if (Root->GraphNr != -1) // if already discovered, step on
1584 Root = Root->next;
1585 }
1586 }
1587 // set cyclic bond criterium on "same LP" basis
1588 Binder = first;
1589 while(Binder->next != last) {
1590 Binder = Binder->next;
1591 if (Binder->rightatom->LowpointNr == Binder->leftatom->LowpointNr) { // cyclic ??
1592 Binder->Cyclic = true;
1593 NoCyclicBonds++;
1594 }
1595 }
1596
1597 // analysis of the cycles (print rings, get minimum cycle length)
1598 CyclicStructureAnalysis(out, BackEdgeStack, MinimumRingSize);
1599
1600 *out << Verbose(1) << "Final graph info for each atom is:" << endl;
1601 Walker = start;
1602 while (Walker->next != end) {
1603 Walker = Walker->next;
1604 *out << Verbose(2) << "Atom " << Walker->Name << " is " << ((Walker->SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
1605 OutputComponentNumber(out, Walker);
1606 *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
1607 }
1608
1609 *out << Verbose(1) << "Final graph info for each bond is:" << endl;
1610 Binder = first;
1611 while(Binder->next != last) {
1612 Binder = Binder->next;
1613 *out << Verbose(2) << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
1614 *out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
1615 OutputComponentNumber(out, Binder->leftatom);
1616 *out << " === ";
1617 *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
1618 OutputComponentNumber(out, Binder->rightatom);
1619 *out << ">." << endl;
1620 if (Binder->Cyclic) // cyclic ??
1621 *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
1622 }
1623
1624 // free all and exit
1625 delete(AtomStack);
1626 *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
1627 return SubGraphs;
1628};
1629
1630/** Analyses the cycles found and returns minimum of all cycle lengths.
1631 * \param *out output stream for debugging
1632 * \param *BackEdgeStack stack with all back edges found during DFS scan
1633 * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
1634 * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
1635 */
1636void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> * BackEdgeStack, int *&MinimumRingSize)
1637{
1638 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CyclicStructureAnalysis: **PredecessorList");
1639 int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
1640 enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
1641 class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount); // will hold the current ring
1642 class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount); // contains all "touched" atoms
1643 atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
1644 bond *Binder = NULL, *BackEdge = NULL;
1645 int RingSize, NumCycles, MinRingSize = -1;
1646
1647 // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
1648 for (int i=0;i<AtomCount;i++) {
1649 PredecessorList[i] = NULL;
1650 ShortestPathList[i] = -1;
1651 ColorList[i] = white;
1652 }
1653 MinimumRingSize = new int[AtomCount];
1654 for(int i=0;i<AtomCount;i++)
1655 MinimumRingSize[i] = AtomCount;
1656
1657
1658 *out << Verbose(1) << "Back edge list - ";
1659 BackEdgeStack->Output(out);
1660
1661 *out << Verbose(1) << "Analysing cycles ... " << endl;
1662 NumCycles = 0;
1663 while (!BackEdgeStack->IsEmpty()) {
1664 BackEdge = BackEdgeStack->PopFirst();
1665 // this is the target
1666 Root = BackEdge->leftatom;
1667 // this is the source point
1668 Walker = BackEdge->rightatom;
1669 ShortestPathList[Walker->nr] = 0;
1670 BFSStack->ClearStack(); // start with empty BFS stack
1671 BFSStack->Push(Walker);
1672 TouchedStack->Push(Walker);
1673 //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
1674 OtherAtom = NULL;
1675 while ((Walker != Root) && ((OtherAtom == NULL) || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]))) { // look for Root
1676 Walker = BFSStack->PopFirst();
1677 //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
1678 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
1679 Binder = ListOfBondsPerAtom[Walker->nr][i];
1680 if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
1681 OtherAtom = Binder->GetOtherAtom(Walker);
1682 //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
1683 if (ColorList[OtherAtom->nr] == white) {
1684 TouchedStack->Push(OtherAtom);
1685 ColorList[OtherAtom->nr] = lightgray;
1686 PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
1687 ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
1688 //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
1689 if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
1690 //*out << Verbose(3) << "Putting OtherAtom into queue." << endl;
1691 BFSStack->Push(OtherAtom);
1692 }
1693 } else {
1694 //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
1695 }
1696 } else {
1697 //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
1698 }
1699 }
1700 ColorList[Walker->nr] = black;
1701 //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
1702 }
1703
1704 if (Walker == Root) {
1705 // now climb back the predecessor list and thus find the cycle members
1706 NumCycles++;
1707 RingSize = 1;
1708 Walker = Root;
1709 *out << Verbose(1) << "Found ring contains: ";
1710 while (Walker != BackEdge->rightatom) {
1711 *out << Walker->Name << " <-> ";
1712 Walker = PredecessorList[Walker->nr];
1713 RingSize++;
1714 }
1715 *out << Walker->Name << " with a length of " << RingSize << "." << endl << endl;
1716 // walk through all and set MinimumRingSize
1717 Walker = Root;
1718 while (Walker != BackEdge->rightatom) {
1719 Walker = PredecessorList[Walker->nr];
1720 if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
1721 MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
1722 }
1723 if ((RingSize < MinRingSize) || (MinRingSize == -1))
1724 MinRingSize = RingSize;
1725 } else {
1726 *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
1727 }
1728
1729 // now clean the lists
1730 while (!TouchedStack->IsEmpty()){
1731 Walker = TouchedStack->PopFirst();
1732 PredecessorList[Walker->nr] = NULL;
1733 ShortestPathList[Walker->nr] = -1;
1734 ColorList[Walker->nr] = white;
1735 }
1736 }
1737 if (MinRingSize != -1) {
1738 // go over all atoms
1739 Root = start;
1740 while(Root->next != end) {
1741 Root = Root->next;
1742
1743 if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
1744 ShortestPathList[Walker->nr] = 0;
1745 BFSStack->ClearStack(); // start with empty BFS stack
1746 BFSStack->Push(Walker);
1747 TouchedStack->Push(Walker);
1748 //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
1749 OtherAtom = Walker;
1750 while ((Walker != Root) && (OtherAtom != NULL)) { // look for Root
1751 Walker = BFSStack->PopFirst();
1752 //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
1753 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
1754 Binder = ListOfBondsPerAtom[Walker->nr][i];
1755 if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
1756 OtherAtom = Binder->GetOtherAtom(Walker);
1757 //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
1758 if (ColorList[OtherAtom->nr] == white) {
1759 TouchedStack->Push(OtherAtom);
1760 ColorList[OtherAtom->nr] = lightgray;
1761 PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
1762 ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
1763 //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
1764 if (MinimumRingSize[OtherAtom->GetTrueFather()->nr] != AtomCount) { // if the other atom is connected to a ring
1765 MinimumRingSize[Root->GetTrueFather()->nr] = ShortestPathList[OtherAtom->nr]+MinimumRingSize[OtherAtom->GetTrueFather()->nr];
1766 OtherAtom = NULL; //break;
1767 break;
1768 } else
1769 BFSStack->Push(OtherAtom);
1770 } else {
1771 //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
1772 }
1773 } else {
1774 //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
1775 }
1776 }
1777 ColorList[Walker->nr] = black;
1778 //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
1779 }
1780
1781 // now clean the lists
1782 while (!TouchedStack->IsEmpty()){
1783 Walker = TouchedStack->PopFirst();
1784 PredecessorList[Walker->nr] = NULL;
1785 ShortestPathList[Walker->nr] = -1;
1786 ColorList[Walker->nr] = white;
1787 }
1788 }
1789 *out << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl;
1790 }
1791 *out << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl;
1792 } else
1793 *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
1794
1795 Free((void **)&PredecessorList, "molecule::CyclicStructureAnalysis: **PredecessorList");
1796 Free((void **)&ShortestPathList, "molecule::CyclicStructureAnalysis: **ShortestPathList");
1797 Free((void **)&ColorList, "molecule::CyclicStructureAnalysis: **ColorList");
1798 delete(BFSStack);
1799};
1800
1801/** Sets the next component number.
1802 * This is O(N) as the number of bonds per atom is bound.
1803 * \param *vertex atom whose next atom::*ComponentNr is to be set
1804 * \param nr number to use
1805 */
1806void molecule::SetNextComponentNumber(atom *vertex, int nr)
1807{
1808 int i=0;
1809 if (vertex != NULL) {
1810 for(;i<NumberOfBondsPerAtom[vertex->nr];i++) {
1811 if (vertex->ComponentNr[i] == -1) { // check if not yet used
1812 vertex->ComponentNr[i] = nr;
1813 break;
1814 }
1815 else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
1816 break; // breaking here will not cause error!
1817 }
1818 if (i == NumberOfBondsPerAtom[vertex->nr])
1819 cerr << "Error: All Component entries are already occupied!" << endl;
1820 } else
1821 cerr << "Error: Given vertex is NULL!" << endl;
1822};
1823
1824/** Output a list of flags, stating whether the bond was visited or not.
1825 * \param *out output stream for debugging
1826 */
1827void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
1828{
1829 for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
1830 *out << vertex->ComponentNr[i] << " ";
1831};
1832
1833/** Allocates memory for all atom::*ComponentNr in this molecule and sets each entry to -1.
1834 */
1835void molecule::InitComponentNumbers()
1836{
1837 atom *Walker = start;
1838 while(Walker->next != end) {
1839 Walker = Walker->next;
1840 if (Walker->ComponentNr != NULL)
1841 Free((void **)&Walker->ComponentNr, "molecule::InitComponentNumbers: **Walker->ComponentNr");
1842 Walker->ComponentNr = (int *) Malloc(sizeof(int)*NumberOfBondsPerAtom[Walker->nr], "molecule::InitComponentNumbers: *Walker->ComponentNr");
1843 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
1844 Walker->ComponentNr[i] = -1;
1845 }
1846};
1847
1848/** Returns next unused bond for this atom \a *vertex or NULL of none exists.
1849 * \param *vertex atom to regard
1850 * \return bond class or NULL
1851 */
1852bond * molecule::FindNextUnused(atom *vertex)
1853{
1854 for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
1855 if (ListOfBondsPerAtom[vertex->nr][i]->IsUsed() == white)
1856 return(ListOfBondsPerAtom[vertex->nr][i]);
1857 return NULL;
1858};
1859
1860/** Resets bond::Used flag of all bonds in this molecule.
1861 * \return true - success, false - -failure
1862 */
1863void molecule::ResetAllBondsToUnused()
1864{
1865 bond *Binder = first;
1866 while (Binder->next != last) {
1867 Binder = Binder->next;
1868 Binder->ResetUsed();
1869 }
1870};
1871
1872/** Resets atom::nr to -1 of all atoms in this molecule.
1873 */
1874void molecule::ResetAllAtomNumbers()
1875{
1876 atom *Walker = start;
1877 while (Walker->next != end) {
1878 Walker = Walker->next;
1879 Walker->GraphNr = -1;
1880 }
1881};
1882
1883/** Output a list of flags, stating whether the bond was visited or not.
1884 * \param *out output stream for debugging
1885 * \param *list
1886 */
1887void OutputAlreadyVisited(ofstream *out, int *list)
1888{
1889 *out << Verbose(4) << "Already Visited Bonds:\t";
1890 for(int i=1;i<=list[0];i++) *out << Verbose(0) << list[i] << " ";
1891 *out << endl;
1892};
1893
1894/** Estimates by educated guessing (using upper limit) the expected number of fragments.
1895 * The upper limit is
1896 * \f[
1897 * n = N \cdot C^k
1898 * \f]
1899 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
1900 * \param *out output stream for debugging
1901 * \param order bond order k
1902 * \return number n of fragments
1903 */
1904int molecule::GuesstimateFragmentCount(ofstream *out, int order)
1905{
1906 int c = 0;
1907 int FragmentCount;
1908 // get maximum bond degree
1909 atom *Walker = start;
1910 while (Walker->next != end) {
1911 Walker = Walker->next;
1912 c = (NumberOfBondsPerAtom[Walker->nr] > c) ? NumberOfBondsPerAtom[Walker->nr] : c;
1913 }
1914 FragmentCount = NoNonHydrogen*(1 << (c*order));
1915 *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
1916 return FragmentCount;
1917};
1918
1919/** Scans a single line for number and puts them into \a KeySet.
1920 * \param *out output stream for debugging
1921 * \param *buffer buffer to scan
1922 * \param &CurrentSet filled KeySet on return
1923 * \return true - at least one valid atom id parsed, false - CurrentSet is empty
1924 */
1925bool molecule::ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
1926{
1927 stringstream line;
1928 int AtomNr;
1929 int status = 0;
1930
1931 line.str(buffer);
1932 while (!line.eof()) {
1933 line >> AtomNr;
1934 if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
1935 CurrentSet.insert(AtomNr); // insert at end, hence in same order as in file!
1936 status++;
1937 } // else it's "-1" or else and thus must not be added
1938 }
1939 *out << Verbose(1) << "The scanned KeySet is ";
1940 for(KeySet::iterator runner = CurrentSet.begin(); runner != CurrentSet.end(); runner++) {
1941 *out << (*runner) << "\t";
1942 }
1943 *out << endl;
1944 return (status != 0);
1945};
1946
1947/** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
1948 * Does two-pass scanning:
1949 * -# Scans the keyset file and initialises a temporary graph
1950 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
1951 * Finally, the temporary graph is inserted into the given \a FragmentList for return.
1952 * \param *out output stream for debugging
1953 * \param *path path to file
1954 * \param *FragmentList empty, filled on return
1955 * \param IsAngstroem whether we have Ansgtroem or bohrradius
1956 * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
1957 */
1958bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList, bool IsAngstroem)
1959{
1960 bool status = true;
1961 ifstream InputFile;
1962 stringstream line;
1963 GraphTestPair testGraphInsert;
1964 int NumberOfFragments = 0;
1965 double TEFactor;
1966 char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
1967
1968 if (FragmentList == NULL) { // check list pointer
1969 FragmentList = new Graph;
1970 }
1971
1972 // 1st pass: open file and read
1973 *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
1974 sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
1975 InputFile.open(filename);
1976 if (InputFile != NULL) {
1977 // each line represents a new fragment
1978 char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - *buffer");
1979 // 1. parse keysets and insert into temp. graph
1980 while (!InputFile.eof()) {
1981 InputFile.getline(buffer, MAXSTRINGSIZE);
1982 KeySet CurrentSet;
1983 if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet))) { // if at least one valid atom was added, write config
1984 testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1))); // store fragment number and current factor
1985 if (!testGraphInsert.second) {
1986 cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
1987 }
1988 //FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
1989 }
1990 }
1991 // 2. Free and done
1992 InputFile.close();
1993 InputFile.clear();
1994 Free((void **)&buffer, "molecule::ParseKeySetFile - *buffer");
1995 *out << Verbose(1) << "done." << endl;
1996 } else {
1997 *out << Verbose(1) << "File " << filename << " not found." << endl;
1998 status = false;
1999 }
2000
2001 // 2nd pass: open TEFactors file and read
2002 *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
2003 sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, TEFACTORSFILE);
2004 InputFile.open(filename);
2005 if (InputFile != NULL) {
2006 // 3. add found TEFactors to each keyset
2007 NumberOfFragments = 0;
2008 for(Graph::iterator runner = FragmentList->begin();runner != FragmentList->end(); runner++) {
2009 if (!InputFile.eof()) {
2010 InputFile >> TEFactor;
2011 (*runner).second.second = TEFactor;
2012 *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
2013 } else {
2014 status = false;
2015 break;
2016 }
2017 }
2018 // 4. Free and done
2019 InputFile.close();
2020 *out << Verbose(1) << "done." << endl;
2021 } else {
2022 *out << Verbose(1) << "File " << filename << " not found." << endl;
2023 status = false;
2024 }
2025
2026 // free memory
2027 Free((void **)&filename, "molecule::ParseKeySetFile - filename");
2028
2029 return status;
2030};
2031
2032/** Stores keysets and TEFactors to file.
2033 * \param *out output stream for debugging
2034 * \param KeySetList Graph with Keysets and factors
2035 * \param *path path to file
2036 * \return true - file written successfully, false - writing failed
2037 */
2038bool molecule::StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
2039{
2040 ofstream output;
2041 bool status = true;
2042 string line;
2043 string::iterator ende;
2044
2045 // open KeySet file
2046 line = path;
2047 line.append("/");
2048 line += FRAGMENTPREFIX;
2049 ende = line.end();
2050 line += KEYSETFILE;
2051 output.open(line.c_str(), ios::out);
2052 *out << Verbose(1) << "Saving key sets of the total graph ... ";
2053 if(output != NULL) {
2054 for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
2055 for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
2056 if (sprinter != (*runner).first.begin())
2057 output << "\t";
2058 output << *sprinter;
2059 }
2060 output << endl;
2061 }
2062 *out << "done." << endl;
2063 } else {
2064 cerr << "Unable to open " << line << " for writing keysets!" << endl;
2065 status = false;
2066 }
2067 output.close();
2068 output.clear();
2069
2070 // open TEFactors file
2071 line.erase(ende, line.end());
2072 line += TEFACTORSFILE;
2073 output.open(line.c_str(), ios::out);
2074 *out << Verbose(1) << "Saving TEFactors of the total graph ... ";
2075 if(output != NULL) {
2076 for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++)
2077 output << (*runner).second.second << endl;
2078 *out << Verbose(1) << "done." << endl;
2079 } else {
2080 *out << Verbose(1) << "failed to open " << line << "." << endl;
2081 status = false;
2082 }
2083 output.close();
2084
2085 return status;
2086};
2087
2088/** Storing the bond structure of a molecule to file.
2089 * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
2090 * \param *out output stream for debugging
2091 * \param *path path to file
2092 * \return true - file written successfully, false - writing failed
2093 */
2094bool molecule::StoreAdjacencyToFile(ofstream *out, char *path)
2095{
2096 ofstream AdjacencyFile;
2097 atom *Walker = NULL;
2098 stringstream line;
2099 bool status = true;
2100
2101 line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
2102 AdjacencyFile.open(line.str().c_str(), ios::out);
2103 *out << Verbose(1) << "Saving adjacency list ... ";
2104 if (AdjacencyFile != NULL) {
2105 Walker = start;
2106 while(Walker->next != end) {
2107 Walker = Walker->next;
2108 AdjacencyFile << Walker->nr << "\t";
2109 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
2110 AdjacencyFile << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
2111 AdjacencyFile << endl;
2112 }
2113 AdjacencyFile.close();
2114 *out << Verbose(1) << "done." << endl;
2115 } else {
2116 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
2117 status = false;
2118 }
2119
2120 return status;
2121};
2122
2123/** Checks contents of adjacency file against bond structure in structure molecule.
2124 * \param *out output stream for debugging
2125 * \param *path path to file
2126 * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
2127 * \return true - structure is equal, false - not equivalence
2128 */
2129bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
2130{
2131 ifstream File;
2132 stringstream line;
2133 bool status = true;
2134 char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
2135
2136 line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
2137 File.open(line.str().c_str(), ios::out);
2138 *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
2139 if (File != NULL) {
2140 // allocate storage structure
2141 int NonMatchNumber = 0; // will number of atoms with differing bond structure
2142 int *CurrentBonds = (int *) Malloc(sizeof(int)*8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
2143 int CurrentBondsOfAtom;
2144
2145 // Parse the file line by line and count the bonds
2146 while (!File.eof()) {
2147 File.getline(buffer, MAXSTRINGSIZE);
2148 stringstream line;
2149 line.str(buffer);
2150 int AtomNr = -1;
2151 line >> AtomNr;
2152 CurrentBondsOfAtom = -1; // we count one too far due to line end
2153 // parse into structure
2154 if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
2155 while (!line.eof())
2156 line >> CurrentBonds[ ++CurrentBondsOfAtom ];
2157 // compare against present bonds
2158 //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
2159 if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
2160 for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
2161 int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
2162 int j = 0;
2163 for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
2164 if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
2165 ListOfAtoms[AtomNr] = NULL;
2166 NonMatchNumber++;
2167 status = false;
2168 //out << "[" << id << "]\t";
2169 } else {
2170 //out << id << "\t";
2171 }
2172 }
2173 //out << endl;
2174 } else {
2175 *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
2176 status = false;
2177 }
2178 }
2179 }
2180 File.close();
2181 File.clear();
2182 if (status) { // if equal we parse the KeySetFile
2183 *out << Verbose(1) << "done: Equal." << endl;
2184 status = true;
2185 } else
2186 *out << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl;
2187 Free((void **)&CurrentBonds, "molecule::CheckAdjacencyFileAgainstMolecule - **CurrentBonds");
2188 } else {
2189 *out << Verbose(1) << "Adjacency file not found." << endl;
2190 status = false;
2191 }
2192 *out << endl;
2193 Free((void **)&buffer, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
2194
2195 return status;
2196};
2197
2198/** Checks whether the OrderAtSite is still bewloe \a Order at some site.
2199 * \param *out output stream for debugging
2200 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
2201 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
2202 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
2203 * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
2204 * \return true - needs further fragmentation, false - does not need fragmentation
2205 */
2206bool molecule::CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, char *path)
2207{
2208 atom *Walker = start;
2209 bool status = false;
2210 ifstream InputFile;
2211
2212 // initialize mask list
2213 for(int i=0;i<AtomCount;i++)
2214 AtomMask[i] = false;
2215
2216 if (Order < 0) { // adaptive increase of BondOrder per site
2217 if (AtomMask[AtomCount] == true) // break after one step
2218 return false;
2219 // parse the EnergyPerFragment file
2220 char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
2221 sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
2222 InputFile.open(buffer, ios::in);
2223 if ((InputFile != NULL) && (GlobalKeySetList != NULL)) {
2224 // transmorph graph keyset list into indexed KeySetList
2225 map<int,KeySet> IndexKeySetList;
2226 for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
2227 IndexKeySetList.insert( pair<int,KeySet>(runner->second.first,runner->first) );
2228 }
2229 int lines = 0;
2230 // count the number of lines, i.e. the number of fragments
2231 InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
2232 InputFile.getline(buffer, MAXSTRINGSIZE);
2233 while(!InputFile.eof()) {
2234 InputFile.getline(buffer, MAXSTRINGSIZE);
2235 lines++;
2236 }
2237 *out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl; // one endline too much
2238 InputFile.clear();
2239 InputFile.seekg(ios::beg);
2240 map<int, pair<double,int> > AdaptiveCriteriaList; // (Root No., (Value, Order)) !
2241 int No, FragOrder;
2242 double Value;
2243 // each line represents a fragment root (Atom::nr) id and its energy contribution
2244 InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
2245 InputFile.getline(buffer, MAXSTRINGSIZE);
2246 while(!InputFile.eof()) {
2247 InputFile.getline(buffer, MAXSTRINGSIZE);
2248 if (strlen(buffer) > 2) {
2249 //*out << Verbose(2) << "Scanning: " << buffer;
2250 stringstream line(buffer);
2251 line >> FragOrder;
2252 line >> ws >> No;
2253 line >> ws >> Value; // skip time entry
2254 line >> ws >> Value;
2255 No -= 1; // indices start at 1 in file, not 0
2256 //*out << Verbose(2) << " - yields (" << No << "," << Value << ")" << endl;
2257
2258 // clean the list of those entries that have been superceded by higher order terms already
2259 map<int,KeySet>::iterator marker = IndexKeySetList.find(No); // find keyset to Frag No.
2260 if (marker != IndexKeySetList.end()) { // if found
2261 // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
2262 pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( Value, Order) ));
2263 map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
2264 if (!InsertedElement.second) { // this root is already present
2265 if ((*PresentItem).second.second < FragOrder) // if order there is lower, update entry with higher-order term
2266 //if ((*PresentItem).second.first < (*runner).first) // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
2267 { // if value is smaller, update value and order
2268 (*PresentItem).second.first = Value;
2269 (*PresentItem).second.second = FragOrder;
2270 *out << Verbose(2) << "Updated element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
2271 }
2272 } else {
2273 *out << Verbose(2) << "Inserted element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
2274 }
2275 } else {
2276 *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
2277 }
2278 }
2279 }
2280 // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
2281 map<double, pair<int,int> > FinalRootCandidates;
2282 *out << Verbose(1) << "Root candidate list is: " << endl;
2283 for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList.begin(); runner != AdaptiveCriteriaList.end(); runner++) {
2284 Walker = FindAtom((*runner).first);
2285 if (Walker != NULL) {
2286 if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
2287 *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
2288 FinalRootCandidates.insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
2289 }
2290 } else {
2291 cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
2292 }
2293 }
2294 // pick the ones still below threshold and mark as to be adaptively updated
2295 for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
2296 No = (*runner).second.first;
2297 *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
2298 AtomMask[No] = true;
2299 status = true;
2300 }
2301 // close and done
2302 InputFile.close();
2303 InputFile.clear();
2304 } else {
2305 cerr << "Unable to parse " << buffer << " file, incrementing all." << endl;
2306 while (Walker->next != end) {
2307 Walker = Walker->next;
2308 #ifdef ADDHYDROGEN
2309 if (Walker->type->Z != 1) // skip hydrogen
2310 #endif
2311 {
2312 AtomMask[Walker->nr] = true; // include all (non-hydrogen) atoms
2313 status = true;
2314 }
2315 }
2316 }
2317 Free((void **)&buffer, "molecule::CheckOrderAtSite: *buffer");
2318 // pick a given number of highest values and set AtomMask
2319 } else { // global increase of Bond Order
2320 while (Walker->next != end) {
2321 Walker = Walker->next;
2322 #ifdef ADDHYDROGEN
2323 if (Walker->type->Z != 1) // skip hydrogen
2324 #endif
2325 {
2326 AtomMask[Walker->nr] = true; // include all (non-hydrogen) atoms
2327 if ((Order != 0) && (Walker->AdaptiveOrder < Order))
2328 status = true;
2329 }
2330 }
2331 if ((Order == 0) && (AtomMask[AtomCount] == true)) // single stepping, just check
2332 status = false;
2333
2334 if (!status)
2335 *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
2336 }
2337
2338 // print atom mask for debugging
2339 *out << " ";
2340 for(int i=0;i<AtomCount;i++)
2341 *out << (i % 10);
2342 *out << endl << "Atom mask is: ";
2343 for(int i=0;i<AtomCount;i++)
2344 *out << (AtomMask[i] ? "t" : "f");
2345 *out << endl;
2346
2347 return status;
2348};
2349
2350/** Create a SortIndex to map from BFS labels to the sequence in which the atoms are given in the config file.
2351 * \param *out output stream for debugging
2352 * \param *&SortIndex Mapping array of size molecule::AtomCount
2353 * \return true - success, false - failure of SortIndex alloc
2354 */
2355bool molecule::CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex)
2356{
2357 element *runner = elemente->start;
2358 int AtomNo = 0;
2359 atom *Walker = NULL;
2360
2361 if (SortIndex != NULL) {
2362 *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
2363 return false;
2364 }
2365 SortIndex = (int *) Malloc(sizeof(int)*AtomCount, "molecule::FragmentMolecule: *SortIndex");
2366 for(int i=0;i<AtomCount;i++)
2367 SortIndex[i] = -1;
2368 while (runner->next != elemente->end) { // go through every element
2369 runner = runner->next;
2370 if (ElementsInMolecule[runner->Z]) { // if this element got atoms
2371 Walker = start;
2372 while (Walker->next != end) { // go through every atom of this element
2373 Walker = Walker->next;
2374 if (Walker->type->Z == runner->Z) // if this atom fits to element
2375 SortIndex[Walker->nr] = AtomNo++;
2376 }
2377 }
2378 }
2379 return true;
2380};
2381
2382/** Performs a many-body bond order analysis for a given bond order.
2383 * -# parses adjacency, keysets and orderatsite files
2384 * -# performs DFS to find connected subgraphs (to leave this in was a design decision: might be useful later)
2385 * -# RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energ
2386y contribution", and that's why this consciously not done in the following loop)
2387 * -# in a loop over all subgraphs
2388 * -# calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
2389 * -# creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)
2390 * -# combines the generated molecule lists from all subgraphs
2391 * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
2392 * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
2393 * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
2394 * subgraph in the MoleculeListClass.
2395 * \param *out output stream for debugging
2396 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
2397 * \param *configuration configuration for writing config files for each fragment
2398 */
2399void molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
2400{
2401 MoleculeListClass *BondFragments = NULL;
2402 int *SortIndex = NULL;
2403 int *MinimumRingSize = NULL;
2404 int FragmentCounter;
2405 MoleculeLeafClass *MolecularWalker = NULL;
2406 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
2407 fstream File;
2408 bool FragmentationToDo = true;
2409 Graph **FragmentList = NULL;
2410 Graph *ParsedFragmentList = NULL;
2411 Graph TotalGraph; // graph with all keysets however local numbers
2412 int TotalNumberOfKeySets = 0;
2413 atom **ListOfAtoms = NULL;
2414 atom ***ListOfLocalAtoms = NULL;
2415 bool *AtomMask = NULL;
2416
2417 *out << endl;
2418#ifdef ADDHYDROGEN
2419 *out << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
2420#else
2421 *out << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl;
2422#endif
2423
2424 // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
2425
2426 // ===== 1. Check whether bond structure is same as stored in files ====
2427
2428 // fill the adjacency list
2429 CreateListOfBondsPerAtom(out);
2430
2431 // create lookup table for Atom::nr
2432 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);
2433
2434 // === compare it with adjacency file ===
2435 FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
2436 Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
2437
2438 // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
2439 Subgraphs = DepthFirstSearchAnalysis(out, false, MinimumRingSize);
2440 // fill the bond structure of the individually stored subgraphs
2441 Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms
2442
2443 // ===== 3. if structure still valid, parse key set file and others =====
2444 FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList, configuration->GetIsAngstroem());
2445
2446 // ===== 4. check globally whether there's something to do actually (first adaptivity check)
2447 FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
2448
2449 // =================================== Begin of FRAGMENTATION ===============================
2450 // ===== 6a. assign each keyset to its respective subgraph =====
2451 Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), false);
2452
2453 KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
2454 AtomMask = new bool[AtomCount+1];
2455 while (CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, configuration->configpath)) {
2456 AtomMask[AtomCount] = true; // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
2457 // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
2458 Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
2459
2460 // ===== 7. fill the bond fragment list =====
2461 FragmentCounter = 0;
2462 MolecularWalker = Subgraphs;
2463 while (MolecularWalker->next != NULL) {
2464 MolecularWalker = MolecularWalker->next;
2465 *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
2466 // output ListOfBondsPerAtom for debugging
2467 MolecularWalker->Leaf->OutputListOfBonds(out);
2468 if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
2469
2470 // call BOSSANOVA method
2471 *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
2472 MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
2473 } else {
2474 cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
2475 }
2476 FragmentCounter++; // next fragment list
2477 }
2478 }
2479 delete[](RootStack);
2480 delete[](AtomMask);
2481 delete(ParsedFragmentList);
2482 delete[](MinimumRingSize);
2483
2484 // free the index lookup list
2485 for (int i=0;i<FragmentCounter;i++)
2486 Free((void **)&ListOfLocalAtoms[i], "molecule::FragmentMolecule - *ListOfLocalAtoms[]");
2487 Free((void **)&ListOfLocalAtoms, "molecule::FragmentMolecule - **ListOfLocalAtoms");
2488
2489 // ==================================== End of FRAGMENTATION ============================================
2490
2491 // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
2492 Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
2493
2494 // free subgraph memory again
2495 FragmentCounter = 0;
2496 if (Subgraphs != NULL) {
2497 while (Subgraphs->next != NULL) {
2498 Subgraphs = Subgraphs->next;
2499 delete(FragmentList[FragmentCounter++]);
2500 delete(Subgraphs->previous);
2501 }
2502 delete(Subgraphs);
2503 }
2504 Free((void **)&FragmentList, "molecule::FragmentMolecule - **FragmentList");
2505
2506 // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
2507 // allocate memory for the pointer array and transmorph graphs into full molecular fragments
2508 BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
2509 int k=0;
2510 for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
2511 KeySet test = (*runner).first;
2512 *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
2513 BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
2514 k++;
2515 }
2516 *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
2517
2518 // ===== 9. Save fragments' configuration and keyset files et al to disk ===
2519 if (BondFragments->NumberOfMolecules != 0) {
2520 // create the SortIndex from BFS labels to order in the config file
2521 CreateMappingLabelsToConfigSequence(out, SortIndex);
2522
2523 *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
2524 if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
2525 *out << Verbose(1) << "All configs written." << endl;
2526 else
2527 *out << Verbose(1) << "Some config writing failed." << endl;
2528
2529 // store force index reference file
2530 BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
2531
2532 // store keysets file
2533 StoreKeySetFile(out, TotalGraph, configuration->configpath);
2534
2535 // store Adjacency file
2536 StoreAdjacencyToFile(out, configuration->configpath);
2537
2538 // store adaptive orders into file
2539 StoreOrderAtSiteFile(out, configuration->configpath);
2540
2541 // restore orbital and Stop values
2542 CalculateOrbitals(*configuration);
2543
2544 // free memory for bond part
2545 *out << Verbose(1) << "Freeing bond memory" << endl;
2546 delete(FragmentList); // remove bond molecule from memory
2547 Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
2548 } else
2549 *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
2550
2551 *out << Verbose(0) << "End of bond fragmentation." << endl;
2552};
2553
2554/** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
2555 * Atoms not present in the file get "-1".
2556 * \param *out output stream for debugging
2557 * \param *path path to file ORDERATSITEFILE
2558 * \return true - file writable, false - not writable
2559 */
2560bool molecule::StoreOrderAtSiteFile(ofstream *out, char *path)
2561{
2562 stringstream line;
2563 ofstream file;
2564
2565 line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
2566 file.open(line.str().c_str());
2567 *out << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
2568 if (file != NULL) {
2569 atom *Walker = start;
2570 while (Walker->next != end) {
2571 Walker = Walker->next;
2572 file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << endl;
2573 *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "." << endl;
2574 }
2575 file.close();
2576 *out << Verbose(1) << "done." << endl;
2577 return true;
2578 } else {
2579 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
2580 return false;
2581 }
2582};
2583
2584/** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
2585 * Atoms not present in the file get "0".
2586 * \param *out output stream for debugging
2587 * \param *path path to file ORDERATSITEFILEe
2588 * \return true - file found and scanned, false - file not found
2589 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
2590 */
2591bool molecule::ParseOrderAtSiteFromFile(ofstream *out, char *path)
2592{
2593 unsigned char *OrderArray = (unsigned char *) Malloc(sizeof(unsigned char)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
2594 bool status;
2595 int AtomNr;
2596 stringstream line;
2597 ifstream file;
2598 int Order;
2599
2600 *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
2601 for(int i=0;i<AtomCount;i++)
2602 OrderArray[i] = 0;
2603 line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
2604 file.open(line.str().c_str());
2605 if (file != NULL) {
2606 for (int i=0;i<AtomCount;i++) // initialise with 0
2607 OrderArray[i] = 0;
2608 while (!file.eof()) { // parse from file
2609 file >> AtomNr;
2610 file >> Order;
2611 OrderArray[AtomNr] = (unsigned char) Order;
2612 //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << "." << endl;
2613 }
2614 atom *Walker = start;
2615 while (Walker->next != end) { // fill into atom classes
2616 Walker = Walker->next;
2617 Walker->AdaptiveOrder = OrderArray[Walker->nr];
2618 *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << "." << endl;
2619 }
2620 file.close();
2621 *out << Verbose(1) << "done." << endl;
2622 status = true;
2623 } else {
2624 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
2625 status = false;
2626 }
2627 Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
2628
2629 *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
2630 return status;
2631};
2632
2633/** Creates an 2d array of pointer with an entry for each atom and each bond it has.
2634 * Updates molecule::ListOfBondsPerAtom, molecule::NumberOfBondsPerAtom by parsing through
2635 * bond chain list, using molecule::AtomCount and molecule::BondCount.
2636 * Allocates memory, fills the array and exits
2637 * \param *out output stream for debugging
2638 */
2639void molecule::CreateListOfBondsPerAtom(ofstream *out)
2640{
2641 bond *Binder = NULL;
2642 atom *Walker = NULL;
2643 int TotalDegree;
2644 *out << Verbose(1) << "Begin of Creating ListOfBondsPerAtom: AtomCount = " << AtomCount << "\tBondCount = " << BondCount << "\tNoNonBonds = " << NoNonBonds << "." << endl;
2645
2646 // re-allocate memory
2647 *out << Verbose(2) << "(Re-)Allocating memory." << endl;
2648 if (ListOfBondsPerAtom != NULL) {
2649 for(int i=0;i<AtomCount;i++)
2650 Free((void **)&ListOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: ListOfBondsPerAtom[i]");
2651 Free((void **)&ListOfBondsPerAtom, "molecule::CreateListOfBondsPerAtom: ListOfBondsPerAtom");
2652 }
2653 if (NumberOfBondsPerAtom != NULL)
2654 Free((void **)&NumberOfBondsPerAtom, "molecule::CreateListOfBondsPerAtom: NumberOfBondsPerAtom");
2655 ListOfBondsPerAtom = (bond ***) Malloc(sizeof(bond **)*AtomCount, "molecule::CreateListOfBondsPerAtom: ***ListOfBondsPerAtom");
2656 NumberOfBondsPerAtom = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfBondsPerAtom: *NumberOfBondsPerAtom");
2657
2658 // reset bond counts per atom
2659 for(int i=0;i<AtomCount;i++)
2660 NumberOfBondsPerAtom[i] = 0;
2661 // count bonds per atom
2662 Binder = first;
2663 while (Binder->next != last) {
2664 Binder = Binder->next;
2665 NumberOfBondsPerAtom[Binder->leftatom->nr]++;
2666 NumberOfBondsPerAtom[Binder->rightatom->nr]++;
2667 }
2668 // allocate list of bonds per atom
2669 for(int i=0;i<AtomCount;i++)
2670 ListOfBondsPerAtom[i] = (bond **) Malloc(sizeof(bond *)*NumberOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: **ListOfBondsPerAtom[]");
2671 // clear the list again, now each NumberOfBondsPerAtom marks current free field
2672 for(int i=0;i<AtomCount;i++)
2673 NumberOfBondsPerAtom[i] = 0;
2674 // fill the list
2675 Binder = first;
2676 while (Binder->next != last) {
2677 Binder = Binder->next;
2678 ListOfBondsPerAtom[Binder->leftatom->nr][NumberOfBondsPerAtom[Binder->leftatom->nr]++] = Binder;
2679 ListOfBondsPerAtom[Binder->rightatom->nr][NumberOfBondsPerAtom[Binder->rightatom->nr]++] = Binder;
2680 }
2681
2682 // output list for debugging
2683 *out << Verbose(3) << "ListOfBondsPerAtom for each atom:" << endl;
2684 Walker = start;
2685 while (Walker->next != end) {
2686 Walker = Walker->next;
2687 *out << Verbose(4) << "Atom " << Walker->Name << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
2688 TotalDegree = 0;
2689 for (int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {
2690 *out << *ListOfBondsPerAtom[Walker->nr][j] << "\t";
2691 TotalDegree += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
2692 }
2693 *out << " -- TotalDegree: " << TotalDegree << endl;
2694 }
2695 *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
2696};
2697
2698/** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
2699 * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
2700 * white and putting into queue.
2701 * \param *out output stream for debugging
2702 * \param *Mol Molecule class to add atoms to
2703 * \param **AddedAtomList list with added atom pointers, index is atom father's number
2704 * \param **AddedBondList list with added bond pointers, index is bond father's number
2705 * \param *Root root vertex for BFS
2706 * \param *Bond bond not to look beyond
2707 * \param BondOrder maximum distance for vertices to add
2708 * \param IsAngstroem lengths are in angstroem or bohrradii
2709 */
2710void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
2711{
2712 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::BreadthFirstSearchAdd: **PredecessorList");
2713 int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::BreadthFirstSearchAdd: *ShortestPathList");
2714 enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::BreadthFirstSearchAdd: *ColorList");
2715 class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
2716 atom *Walker = NULL, *OtherAtom = NULL;
2717 bond *Binder = NULL;
2718
2719 // add Root if not done yet
2720 AtomStack->ClearStack();
2721 if (AddedAtomList[Root->nr] == NULL) // add Root if not yet present
2722 AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
2723 AtomStack->Push(Root);
2724
2725 // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
2726 for (int i=0;i<AtomCount;i++) {
2727 PredecessorList[i] = NULL;
2728 ShortestPathList[i] = -1;
2729 if (AddedAtomList[i] != NULL) // mark already present atoms (i.e. Root and maybe others) as visited
2730 ColorList[i] = lightgray;
2731 else
2732 ColorList[i] = white;
2733 }
2734 ShortestPathList[Root->nr] = 0;
2735
2736 // and go on ... Queue always contains all lightgray vertices
2737 while (!AtomStack->IsEmpty()) {
2738 // we have to pop the oldest atom from stack. This keeps the atoms on the stack always of the same ShortestPath distance.
2739 // e.g. if current atom is 2, push to end of stack are of length 3, but first all of length 2 would be popped. They again
2740 // append length of 3 (their neighbours). Thus on stack we have always atoms of a certain length n at bottom of stack and
2741 // followed by n+1 till top of stack.
2742 Walker = AtomStack->PopFirst(); // pop oldest added
2743 *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
2744 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
2745 Binder = ListOfBondsPerAtom[Walker->nr][i];
2746 if (Binder != NULL) { // don't look at bond equal NULL
2747 OtherAtom = Binder->GetOtherAtom(Walker);
2748 *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
2749 if (ColorList[OtherAtom->nr] == white) {
2750 if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
2751 ColorList[OtherAtom->nr] = lightgray;
2752 PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
2753 ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
2754 *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
2755 if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
2756 *out << Verbose(3);
2757 if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
2758 AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
2759 *out << "Added OtherAtom " << OtherAtom->Name;
2760 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2761 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2762 AddedBondList[Binder->nr]->Type = Binder->Type;
2763 *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
2764 } else { // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
2765 *out << "Not adding OtherAtom " << OtherAtom->Name;
2766 if (AddedBondList[Binder->nr] == NULL) {
2767 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2768 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2769 AddedBondList[Binder->nr]->Type = Binder->Type;
2770 *out << ", added Bond " << *(AddedBondList[Binder->nr]);
2771 } else
2772 *out << ", not added Bond ";
2773 }
2774 *out << ", putting OtherAtom into queue." << endl;
2775 AtomStack->Push(OtherAtom);
2776 } else { // out of bond order, then replace
2777 if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
2778 ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
2779 if (Binder == Bond)
2780 *out << Verbose(3) << "Not Queueing, is the Root bond";
2781 else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
2782 *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
2783 if (!Binder->Cyclic)
2784 *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
2785 if (AddedBondList[Binder->nr] == NULL) {
2786 if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
2787 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2788 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2789 AddedBondList[Binder->nr]->Type = Binder->Type;
2790 } else {
2791#ifdef ADDHYDROGEN
2792 Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem);
2793#endif
2794 }
2795 }
2796 }
2797 } else {
2798 *out << Verbose(3) << "Not Adding, has already been visited." << endl;
2799 // This has to be a cyclic bond, check whether it's present ...
2800 if (AddedBondList[Binder->nr] == NULL) {
2801 if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
2802 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2803 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2804 AddedBondList[Binder->nr]->Type = Binder->Type;
2805 } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
2806#ifdef ADDHYDROGEN
2807 Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem);
2808#endif
2809 }
2810 }
2811 }
2812 }
2813 }
2814 ColorList[Walker->nr] = black;
2815 *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
2816 }
2817 Free((void **)&PredecessorList, "molecule::BreadthFirstSearchAdd: **PredecessorList");
2818 Free((void **)&ShortestPathList, "molecule::BreadthFirstSearchAdd: **ShortestPathList");
2819 Free((void **)&ColorList, "molecule::BreadthFirstSearchAdd: **ColorList");
2820 delete(AtomStack);
2821};
2822
2823/** Adds bond structure to this molecule from \a Father molecule.
2824 * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
2825 * with end points present in this molecule, bond is created in this molecule.
2826 * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
2827 * \param *out output stream for debugging
2828 * \param *Father father molecule
2829 * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
2830 * \todo not checked, not fully working probably
2831 */
2832bool molecule::BuildInducedSubgraph(ofstream *out, const molecule *Father)
2833{
2834 atom *Walker = NULL, *OtherAtom = NULL;
2835 bool status = true;
2836 atom **ParentList = (atom **) Malloc(sizeof(atom *)*Father->AtomCount, "molecule::BuildInducedSubgraph: **ParentList");
2837
2838 *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
2839
2840 // reset parent list
2841 *out << Verbose(3) << "Resetting ParentList." << endl;
2842 for (int i=0;i<Father->AtomCount;i++)
2843 ParentList[i] = NULL;
2844
2845 // fill parent list with sons
2846 *out << Verbose(3) << "Filling Parent List." << endl;
2847 Walker = start;
2848 while (Walker->next != end) {
2849 Walker = Walker->next;
2850 ParentList[Walker->father->nr] = Walker;
2851 // Outputting List for debugging
2852 *out << Verbose(4) << "Son["<< Walker->father->nr <<"] of " << Walker->father << " is " << ParentList[Walker->father->nr] << "." << endl;
2853 }
2854
2855 // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
2856 *out << Verbose(3) << "Creating bonds." << endl;
2857 Walker = Father->start;
2858 while (Walker->next != Father->end) {
2859 Walker = Walker->next;
2860 if (ParentList[Walker->nr] != NULL) {
2861 if (ParentList[Walker->nr]->father != Walker) {
2862 status = false;
2863 } else {
2864 for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
2865 OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
2866 if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
2867 *out << Verbose(4) << "Endpoints of Bond " << Father->ListOfBondsPerAtom[Walker->nr][i] << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
2868 AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], Father->ListOfBondsPerAtom[Walker->nr][i]->BondDegree);
2869 }
2870 }
2871 }
2872 }
2873 }
2874
2875 Free((void **)&ParentList, "molecule::BuildInducedSubgraph: **ParentList");
2876 *out << Verbose(2) << "End of BuildInducedSubgraph." << endl;
2877 return status;
2878};
2879
2880
2881/** Looks through a StackClass<atom *> and returns the likeliest removal candiate.
2882 * \param *out output stream for debugging messages
2883 * \param *&Leaf KeySet to look through
2884 * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
2885 * \param index of the atom suggested for removal
2886 */
2887int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
2888{
2889 atom *Runner = NULL;
2890 int SP, Removal;
2891
2892 *out << Verbose(2) << "Looking for removal candidate." << endl;
2893 SP = -1; //0; // not -1, so that Root is never removed
2894 Removal = -1;
2895 for (KeySet::iterator runner = Leaf->begin(); runner != Leaf->end(); runner++) {
2896 Runner = FindAtom((*runner));
2897 if (Runner->type->Z != 1) { // skip all those added hydrogens when re-filling snake stack
2898 if (ShortestPathList[(*runner)] > SP) { // remove the oldest one with longest shortest path
2899 SP = ShortestPathList[(*runner)];
2900 Removal = (*runner);
2901 }
2902 }
2903 }
2904 return Removal;
2905};
2906
2907/** Stores a fragment from \a KeySet into \a molecule.
2908 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
2909 * molecule and adds missing hydrogen where bonds were cut.
2910 * \param *out output stream for debugging messages
2911 * \param &Leaflet pointer to KeySet structure
2912 * \param IsAngstroem whether we have Ansgtroem or bohrradius
2913 * \return pointer to constructed molecule
2914 */
2915molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
2916{
2917 atom *Runner = NULL, *FatherOfRunner = NULL, *OtherFather = NULL;
2918 atom **SonList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::StoreFragmentFromStack: **SonList");
2919 molecule *Leaf = new molecule(elemente);
2920
2921// *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
2922
2923 Leaf->BondDistance = BondDistance;
2924 for(int i=0;i<NDIM*2;i++)
2925 Leaf->cell_size[i] = cell_size[i];
2926
2927 // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
2928 for(int i=0;i<AtomCount;i++)
2929 SonList[i] = NULL;
2930
2931 // first create the minimal set of atoms from the KeySet
2932 for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
2933 FatherOfRunner = FindAtom((*runner)); // find the id
2934 SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
2935 }
2936
2937 // create the bonds between all: Make it an induced subgraph and add hydrogen
2938// *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
2939 Runner = Leaf->start;
2940 while (Runner->next != Leaf->end) {
2941 Runner = Runner->next;
2942 FatherOfRunner = Runner->father;
2943 if (SonList[FatherOfRunner->nr] != NULL) { // check if this, our father, is present in list
2944 // create all bonds
2945 for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
2946 OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
2947// *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
2948 if (SonList[OtherFather->nr] != NULL) {
2949// *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
2950 if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
2951// *out << Verbose(3) << "Adding Bond: " << Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree) << "." << endl;
2952 //NumBonds[Runner->nr]++;
2953 } else {
2954// *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
2955 }
2956 } else {
2957// *out << ", who has no son in this fragment molecule." << endl;
2958#ifdef ADDHYDROGEN
2959// *out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
2960 Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem);
2961#endif
2962 //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
2963 }
2964 }
2965 } else {
2966 *out << Verbose(0) << "ERROR: Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
2967 }
2968#ifdef ADDHYDROGEN
2969 while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
2970 Runner = Runner->next;
2971#endif
2972 }
2973 Leaf->CreateListOfBondsPerAtom(out);
2974 //Leaflet->Leaf->ScanForPeriodicCorrection(out);
2975 Free((void **)&SonList, "molecule::StoreFragmentFromStack: **SonList");
2976// *out << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
2977 return Leaf;
2978};
2979
2980/** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
2981 * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
2982 * computer game, that winds through the connected graph representing the molecule. Color (white,
2983 * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
2984 * creating only unique fragments and not additional ones with vertices simply in different sequence.
2985 * The Predecessor is always the one that came before in discovering, needed on backstepping. And
2986 * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
2987 * stepping.
2988 * \param *out output stream for debugging
2989 * \param Order number of atoms in each fragment
2990 * \param *configuration configuration for writing config files for each fragment
2991 * \return List of all unique fragments with \a Order atoms
2992 */
2993/*
2994MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
2995{
2996 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
2997 int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
2998 int *Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
2999 enum Shading *ColorVertexList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
3000 enum Shading *ColorEdgeList = (enum Shading *) Malloc(sizeof(enum Shading)*BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
3001 StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
3002 StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
3003 StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
3004 MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
3005 MoleculeListClass *FragmentList = NULL;
3006 atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
3007 bond *Binder = NULL;
3008 int RunningIndex = 0, FragmentCounter = 0;
3009
3010 *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
3011
3012 // reset parent list
3013 *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
3014 for (int i=0;i<AtomCount;i++) { // reset all atom labels
3015 // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
3016 Labels[i] = -1;
3017 SonList[i] = NULL;
3018 PredecessorList[i] = NULL;
3019 ColorVertexList[i] = white;
3020 ShortestPathList[i] = -1;
3021 }
3022 for (int i=0;i<BondCount;i++)
3023 ColorEdgeList[i] = white;
3024 RootStack->ClearStack(); // clearstack and push first atom if exists
3025 TouchedStack->ClearStack();
3026 Walker = start->next;
3027 while ((Walker != end)
3028#ifdef ADDHYDROGEN
3029 && (Walker->type->Z == 1)
3030#endif
3031 ) { // search for first non-hydrogen atom
3032 *out << Verbose(4) << "Current Root candidate is " << Walker->Name << "." << endl;
3033 Walker = Walker->next;
3034 }
3035 if (Walker != end)
3036 RootStack->Push(Walker);
3037 else
3038 *out << Verbose(0) << "ERROR: Could not find an appropriate Root atom!" << endl;
3039 *out << Verbose(3) << "Root " << Walker->Name << " is on AtomStack, beginning loop through all vertices ..." << endl;
3040
3041 ///// OUTER LOOP ////////////
3042 while (!RootStack->IsEmpty()) {
3043 // get new root vertex from atom stack
3044 Root = RootStack->PopFirst();
3045 ShortestPathList[Root->nr] = 0;
3046 if (Labels[Root->nr] == -1)
3047 Labels[Root->nr] = RunningIndex++; // prevent it from getting again on AtomStack
3048 PredecessorList[Root->nr] = Root;
3049 TouchedStack->Push(Root);
3050 *out << Verbose(0) << "Root for this loop is: " << Root->Name << ".\n";
3051
3052 // clear snake stack
3053 SnakeStack->ClearStack();
3054 //SnakeStack->TestImplementation(out, start->next);
3055
3056 ///// INNER LOOP ////////////
3057 // Problems:
3058 // - what about cyclic bonds?
3059 Walker = Root;
3060 do {
3061 *out << Verbose(1) << "Current Walker is: " << Walker->Name;
3062 // initial setting of the new Walker: label, color, shortest path and put on stacks
3063 if (Labels[Walker->nr] == -1) { // give atom a unique, monotonely increasing number
3064 Labels[Walker->nr] = RunningIndex++;
3065 RootStack->Push(Walker);
3066 }
3067 *out << ", has label " << Labels[Walker->nr];
3068 if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) { // color it if newly discovered and push on stacks (and if within reach!)
3069 if ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white)) {
3070 // Binder ought to be set still from last neighbour search
3071 *out << ", coloring bond " << *Binder << " black";
3072 ColorEdgeList[Binder->nr] = black; // mark this bond as used
3073 }
3074 if (ShortestPathList[Walker->nr] == -1) {
3075 ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1;
3076 TouchedStack->Push(Walker); // mark every atom for lists cleanup later, whose shortest path has been changed
3077 }
3078 if ((ShortestPathList[Walker->nr] < Order) && (ColorVertexList[Walker->nr] != darkgray)) { // if not already on snake stack
3079 SnakeStack->Push(Walker);
3080 ColorVertexList[Walker->nr] = darkgray; // mark as dark gray of on snake stack
3081 }
3082 }
3083 *out << ", SP of " << ShortestPathList[Walker->nr] << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
3084
3085 // then check the stack for a newly stumbled upon fragment
3086 if (SnakeStack->ItemCount() == Order) { // is stack full?
3087 // store the fragment if it is one and get a removal candidate
3088 Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
3089 // remove the candidate if one was found
3090 if (Removal != NULL) {
3091 *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
3092 SnakeStack->RemoveItem(Removal);
3093 ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
3094 if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
3095 Walker = PredecessorList[Removal->nr];
3096 *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
3097 }
3098 }
3099 } else
3100 Removal = NULL;
3101
3102 // finally, look for a white neighbour as the next Walker
3103 Binder = NULL;
3104 if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) { // don't look, if a new walker has been set above
3105 *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
3106 OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
3107 if (ShortestPathList[Walker->nr] < Order) {
3108 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
3109 Binder = ListOfBondsPerAtom[Walker->nr][i];
3110 *out << Verbose(2) << "Current bond is " << *Binder << ": ";
3111 OtherAtom = Binder->GetOtherAtom(Walker);
3112 if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
3113 *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
3114 //ColorVertexList[OtherAtom->nr] = lightgray; // mark as explored
3115 } else { // otherwise check its colour and element
3116 if (
3117#ifdef ADDHYDROGEN
3118 (OtherAtom->type->Z != 1) &&
3119#endif
3120 (ColorEdgeList[Binder->nr] == white)) { // skip hydrogen, look for unexplored vertices
3121 *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
3122 // i find it currently rather sensible to always set the predecessor in order to find one's way back
3123 //if (PredecessorList[OtherAtom->nr] == NULL) {
3124 PredecessorList[OtherAtom->nr] = Walker;
3125 *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
3126 //} else {
3127 // *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
3128 //}
3129 Walker = OtherAtom;
3130 break;
3131 } else {
3132 if (OtherAtom->type->Z == 1)
3133 *out << "Links to a hydrogen atom." << endl;
3134 else
3135 *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
3136 }
3137 }
3138 }
3139 } else { // means we have stepped beyond the horizon: Return!
3140 Walker = PredecessorList[Walker->nr];
3141 OtherAtom = Walker;
3142 *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
3143 }
3144 if (Walker != OtherAtom) { // if no white neighbours anymore, color it black
3145 *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
3146 ColorVertexList[Walker->nr] = black;
3147 Walker = PredecessorList[Walker->nr];
3148 }
3149 }
3150 } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
3151 *out << Verbose(2) << "Inner Looping is finished." << endl;
3152
3153 // if we reset all AtomCount atoms, we have again technically O(N^2) ...
3154 *out << Verbose(2) << "Resetting lists." << endl;
3155 Walker = NULL;
3156 Binder = NULL;
3157 while (!TouchedStack->IsEmpty()) {
3158 Walker = TouchedStack->PopLast();
3159 *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
3160 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
3161 ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
3162 PredecessorList[Walker->nr] = NULL;
3163 ColorVertexList[Walker->nr] = white;
3164 ShortestPathList[Walker->nr] = -1;
3165 }
3166 }
3167 *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
3168
3169 // copy together
3170 *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
3171 FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
3172 RunningIndex = 0;
3173 while ((Leaflet != NULL) && (RunningIndex < FragmentCounter)) {
3174 FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
3175 Leaflet->Leaf = NULL; // prevent molecule from being removed
3176 TempLeaf = Leaflet;
3177 Leaflet = Leaflet->previous;
3178 delete(TempLeaf);
3179 };
3180
3181 // free memory and exit
3182 Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
3183 Free((void **)&ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
3184 Free((void **)&Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
3185 Free((void **)&ColorVertexList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
3186 delete(RootStack);
3187 delete(TouchedStack);
3188 delete(SnakeStack);
3189
3190 *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
3191 return FragmentList;
3192};
3193*/
3194
3195/** Structure containing all values in power set combination generation.
3196 */
3197struct UniqueFragments {
3198 config *configuration;
3199 atom *Root;
3200 Graph *Leaflet;
3201 KeySet *FragmentSet;
3202 int ANOVAOrder;
3203 int FragmentCounter;
3204 int CurrentIndex;
3205 int *Labels;
3206 double TEFactor;
3207 int *ShortestPathList;
3208 bool **UsedList;
3209 bond **BondsPerSPList;
3210 int *BondsPerSPCount;
3211};
3212
3213/** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.
3214 * -# loops over every possible combination (2^dimension of edge set)
3215 * -# inserts current set, if there's still space left
3216 * -# yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root dist
3217ance+1
3218 * -# no: stores fragment into keyset list by calling InsertFragmentIntoGraph
3219 * -# removes all items added into the snake stack (in UniqueFragments structure) added during level (root
3220distance) and current set
3221 * \param *out output stream for debugging
3222 * \param FragmentSearch UniqueFragments structure with all values needed
3223 * \param RootDistance current shortest path level, whose set of edges is represented by **BondsSet
3224 * \param SetDimension Number of possible bonds on this level (i.e. size of the array BondsSet[])
3225 * \param SubOrder remaining number of allowed vertices to add
3226 */
3227void molecule::SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder)
3228{
3229 atom *OtherWalker = NULL;
3230 int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
3231 int NumCombinations;
3232 bool bit;
3233 int bits, TouchedIndex, SubSetDimension, SP;
3234 int Removal;
3235 int *TouchedList = (int *) Malloc(sizeof(int)*(SubOrder+1), "molecule::SPFragmentGenerator: *TouchedList");
3236 bond *Binder = NULL;
3237 bond **BondsList = NULL;
3238
3239 NumCombinations = 1 << SetDimension;
3240
3241 // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
3242 // von Endstuecken (aus den Bonds) hinzugefÃŒgt werden und fÃŒr verbleibende ANOVAOrder
3243 // rekursiv GraphCrawler in der nÀchsten Ebene aufgerufen werden
3244
3245 *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
3246 *out << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->Root << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " << NumCombinations-1 << " combination(s)." << endl;
3247
3248 // initialised touched list (stores added atoms on this level)
3249 *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
3250 for (TouchedIndex=0;TouchedIndex<=SubOrder;TouchedIndex++) // empty touched list
3251 TouchedList[TouchedIndex] = -1;
3252 TouchedIndex = 0;
3253
3254 // create every possible combination of the endpieces
3255 *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
3256 for (int i=1;i<NumCombinations;i++) { // sweep through all power set combinations (skip empty set!)
3257 // count the set bit of i
3258 bits = 0;
3259 for (int j=0;j<SetDimension;j++)
3260 bits += (i & (1 << j)) >> j;
3261
3262 *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
3263 if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
3264 // --1-- add this set of the power set of bond partners to the snake stack
3265 for (int j=0;j<SetDimension;j++) { // pull out every bit by shifting
3266 bit = ((i & (1 << j)) != 0); // mask the bit for the j-th bond
3267 if (bit) { // if bit is set, we add this bond partner
3268 OtherWalker = BondsSet[j]->rightatom; // rightatom is always the one more distant, i.e. the one to add
3269 //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
3270 //if ((!FragmentSearch->UsedList[OtherWalker->nr][i]) && (FragmentSearch->Labels[OtherWalker->nr] > FragmentSearch->Labels[FragmentSearch->Root->nr])) {
3271 //*out << Verbose(2+verbosity) << "Not used so far, label " << FragmentSearch->Labels[OtherWalker->nr] << " is bigger than Root's " << FragmentSearch->Labels[FragmentSearch->Root->nr] << "." << endl;
3272 *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
3273 TouchedList[TouchedIndex++] = OtherWalker->nr; // note as added
3274 FragmentSearch->FragmentSet->insert(OtherWalker->nr);
3275 //FragmentSearch->UsedList[OtherWalker->nr][i] = true;
3276 //}
3277 } else {
3278 *out << Verbose(2+verbosity) << "Not adding." << endl;
3279 }
3280 }
3281
3282 if (bits < SubOrder) {
3283 *out << Verbose(1+verbosity) << "There's still some space left on stack: " << (SubOrder - bits) << "." << endl;
3284 // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
3285 SP = RootDistance+1; // this is the next level
3286 // first count the members in the subset
3287 SubSetDimension = 0;
3288 Binder = FragmentSearch->BondsPerSPList[2*SP]; // start node for this level
3289 while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) { // compare to end node of this level
3290 Binder = Binder->next;
3291 for (int k=0;k<TouchedIndex;k++) {
3292 if (Binder->Contains(TouchedList[k])) // if we added this very endpiece
3293 SubSetDimension++;
3294 }
3295 }
3296 // then allocate and fill the list
3297 BondsList = (bond **) Malloc(sizeof(bond *)*SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
3298 SubSetDimension = 0;
3299 Binder = FragmentSearch->BondsPerSPList[2*SP];
3300 while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {
3301 Binder = Binder->next;
3302 for (int k=0;k<TouchedIndex;k++) {
3303 if (Binder->leftatom->nr == TouchedList[k]) // leftatom is always the close one
3304 BondsList[SubSetDimension++] = Binder;
3305 }
3306 }
3307 *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
3308 SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
3309 Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
3310 } else {
3311 // --2-- otherwise store the complete fragment
3312 *out << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
3313 // store fragment as a KeySet
3314 *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
3315 for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
3316 *out << (*runner) << " ";
3317 InsertFragmentIntoGraph(out, FragmentSearch);
3318 //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
3319 //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList,FragmentSearch->Labels, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
3320 }
3321
3322 // --3-- remove all added items in this level from snake stack
3323 *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
3324 for(int j=0;j<TouchedIndex;j++) {
3325 Removal = TouchedList[j];
3326 *out << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
3327 FragmentSearch->FragmentSet->erase(Removal);
3328 TouchedList[j] = -1;
3329 }
3330 *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
3331 for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
3332 *out << (*runner) << " ";
3333 *out << endl;
3334 TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
3335 } else {
3336 *out << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
3337 }
3338 }
3339 Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");
3340 *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
3341};
3342
3343/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
3344 * -# initialises UniqueFragments structure
3345 * -# fills edge list via BFS
3346 * -# creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as
3347 root distance, the edge set, its dimension and the current suborder
3348 * -# Free'ing structure
3349 * Note that we may use the fact that the atoms are SP-ordered on the atomstack. I.e. when popping always the last, we first get all
3350 * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
3351 * \param *out output stream for debugging
3352 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
3353 * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
3354 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
3355 * \return number of inserted fragments
3356 * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
3357 */
3358int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
3359{
3360 int SP, UniqueIndex, AtomKeyNr;
3361 int *NumberOfAtomsSPLevel = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *SPLevelCount");
3362 atom *Walker = NULL, *OtherWalker = NULL;
3363 bond *Binder = NULL;
3364 bond **BondsList = NULL;
3365 KeyStack AtomStack;
3366 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
3367 KeySet::iterator runner;
3368 int RootKeyNr = FragmentSearch.Root->nr;
3369 int Counter = FragmentSearch.FragmentCounter;
3370
3371 for (int i=0;i<AtomCount;i++) {
3372 PredecessorList[i] = NULL;
3373 }
3374
3375 *out << endl;
3376 *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;
3377
3378 UniqueIndex = 0;
3379 if (FragmentSearch.Labels[RootKeyNr] == -1)
3380 FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
3381 FragmentSearch.ShortestPathList[RootKeyNr] = 0;
3382 // prepare the atom stack counters (number of atoms with certain SP on stack)
3383 for (int i=0;i<Order;i++)
3384 NumberOfAtomsSPLevel[i] = 0;
3385 NumberOfAtomsSPLevel[0] = 1; // for root
3386 SP = -1;
3387 *out << endl;
3388 *out << Verbose(0) << "Starting BFS analysis ..." << endl;
3389 // push as first on atom stack and goooo ...
3390 AtomStack.clear();
3391 AtomStack.push_back(RootKeyNr);
3392 *out << Verbose(0) << "Cleared AtomStack and pushed root as first item onto it." << endl;
3393 // do a BFS search to fill the SP lists and label the found vertices
3394 while (!AtomStack.empty()) {
3395 // pop next atom
3396 AtomKeyNr = AtomStack.front();
3397 AtomStack.pop_front();
3398 if (SP != -1)
3399 NumberOfAtomsSPLevel[SP]--;
3400 if ((SP == -1) || (NumberOfAtomsSPLevel[SP] == -1)) {
3401 SP++;
3402 NumberOfAtomsSPLevel[SP]--; // carry over "-1" to next level
3403 *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with 0 item(s)";
3404 if (SP > 0)
3405 *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
3406 else
3407 *out << "." << endl;
3408 FragmentSearch.BondsPerSPCount[SP] = 0;
3409 } else {
3410 *out << Verbose(1) << "Still " << NumberOfAtomsSPLevel[SP]+1 << " on this SP (" << SP << ") level, currently having " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
3411 }
3412 Walker = FindAtom(AtomKeyNr);
3413 *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << "." << endl;
3414 // check for new sp level
3415 // go through all its bonds
3416 *out << Verbose(1) << "Going through all bonds of Walker." << endl;
3417 for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
3418 Binder = ListOfBondsPerAtom[AtomKeyNr][i];
3419 OtherWalker = Binder->GetOtherAtom(Walker);
3420 if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
3421#ifdef ADDHYDROGEN
3422 && (OtherWalker->type->Z != 1)
3423#endif
3424 ) { // skip hydrogens and restrict to fragment
3425 *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
3426 // set the label if not set (and push on root stack as well)
3427 if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
3428 FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
3429 *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
3430 } else {
3431 *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
3432 }
3433 if ((OtherWalker != PredecessorList[AtomKeyNr]) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
3434 FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
3435 *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
3436 if (FragmentSearch.ShortestPathList[OtherWalker->nr] < Order) { // only pass through those within reach of Order
3437 // push for exploration on stack (only if the SP of OtherWalker is longer than of Walker! (otherwise it has been added already!)
3438 if (FragmentSearch.ShortestPathList[OtherWalker->nr] > SP) {
3439 if (FragmentSearch.ShortestPathList[OtherWalker->nr] < (Order-1)) {
3440 *out << Verbose(3) << "Putting on atom stack for further exploration." << endl;
3441 PredecessorList[OtherWalker->nr] = Walker; // note down the one further up the exploration tree
3442 AtomStack.push_back(OtherWalker->nr);
3443 NumberOfAtomsSPLevel[FragmentSearch.ShortestPathList[OtherWalker->nr]]++;
3444 } else {
3445 *out << Verbose(3) << "Not putting on atom stack, is at end of reach." << endl;
3446 }
3447 // add the bond in between to the SP list
3448 Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
3449 add(Binder, FragmentSearch.BondsPerSPList[2*SP+1]);
3450 FragmentSearch.BondsPerSPCount[SP]++;
3451 *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
3452 } else *out << Verbose(3) << "Not putting on atom stack, nor adding bond, as " << *OtherWalker << "'s SP is shorter than Walker's." << endl;
3453 } else *out << Verbose(3) << "Not continuing, as " << *OtherWalker << " is out of reach of order " << Order << "." << endl;
3454 } else *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->nr << " is smaller than that of Root " << RootKeyNr << " or this is my predecessor." << endl;
3455 } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
3456 }
3457 }
3458 // reset predecessor list
3459 for(int i=0;i<Order;i++) {
3460 Binder = FragmentSearch.BondsPerSPList[2*i];
3461 *out << Verbose(1) << "Current SP level is " << i << "." << endl;
3462 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3463 Binder = Binder->next;
3464 PredecessorList[Binder->rightatom->nr] = NULL; // By construction "OtherAtom" is always Bond::rightatom!
3465 }
3466 }
3467 *out << endl;
3468
3469 // outputting all list for debugging
3470 *out << Verbose(0) << "Printing all found lists." << endl;
3471 for(int i=0;i<Order;i++) {
3472 Binder = FragmentSearch.BondsPerSPList[2*i];
3473 *out << Verbose(1) << "Current SP level is " << i << "." << endl;
3474 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3475 Binder = Binder->next;
3476 *out << Verbose(2) << *Binder << endl;
3477 }
3478 }
3479
3480 // creating fragments with the found edge sets
3481 SP = 0;
3482 for(int i=0;i<Order;i++) { // sum up all found edges
3483 Binder = FragmentSearch.BondsPerSPList[2*i];
3484 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3485 Binder = Binder->next;
3486 SP ++;
3487 }
3488 }
3489 *out << Verbose(0) << "Total number of edges is " << SP << "." << endl;
3490 if (SP >= (Order-1)) {
3491 // start with root (push on fragment stack)
3492 *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
3493 FragmentSearch.FragmentSet->clear();
3494 FragmentSearch.FragmentSet->insert(FragmentSearch.Root->nr);
3495
3496 if (FragmentSearch.FragmentSet->size() == (unsigned int) Order) {
3497 *out << Verbose(0) << "Enough items on stack already for a fragment!" << endl;
3498 // store fragment as a KeySet
3499 *out << Verbose(2) << "Found a new fragment[" << FragmentSearch.FragmentCounter << "], local nr.s are: ";
3500 for(KeySet::iterator runner = FragmentSearch.FragmentSet->begin(); runner != FragmentSearch.FragmentSet->end(); runner++) {
3501 *out << (*runner) << " ";
3502 }
3503 *out << endl;
3504 InsertFragmentIntoGraph(out, &FragmentSearch);
3505 } else {
3506 *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
3507 // prepare the subset and call the generator
3508 BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
3509 Binder = FragmentSearch.BondsPerSPList[0];
3510 for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++) {
3511 Binder = Binder->next;
3512 BondsList[i] = Binder;
3513 }
3514 SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order-1);
3515 Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
3516 }
3517 } else {
3518 *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
3519 }
3520
3521 // as FragmentSearch structure is used only once, we don't have to clean it anymore
3522 // remove root from stack
3523 *out << Verbose(0) << "Removing root again from stack." << endl;
3524 FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
3525
3526 // free'ing the bonds lists
3527 *out << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl;
3528 for(int i=0;i<Order;i++) {
3529 *out << Verbose(1) << "Current SP level is " << i << ": ";
3530 Binder = FragmentSearch.BondsPerSPList[2*i];
3531 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3532 Binder = Binder->next;
3533 // *out << "Removing atom " << Binder->leftatom->nr << " and " << Binder->rightatom->nr << "." << endl; // make sure numbers are local
3534 FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
3535 FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
3536 }
3537 // delete added bonds
3538 cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
3539 // also start and end node
3540 *out << "cleaned." << endl;
3541 }
3542
3543 // free allocated memory
3544 Free((void **)&NumberOfAtomsSPLevel, "molecule::PowerSetGenerator: *SPLevelCount");
3545 Free((void **)&PredecessorList, "molecule::PowerSetGenerator: **PredecessorList");
3546
3547 // return list
3548 *out << Verbose(0) << "End of PowerSetGenerator." << endl;
3549 return (FragmentSearch.FragmentCounter - Counter);
3550};
3551
3552/** Corrects the nuclei position if the fragment was created over the cell borders.
3553 * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
3554 * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
3555 * and re-add the bond. Looping on the distance check.
3556 * \param *out ofstream for debugging messages
3557 */
3558void molecule::ScanForPeriodicCorrection(ofstream *out)
3559{
3560 bond *Binder = NULL;
3561 bond *OtherBinder = NULL;
3562 atom *Walker = NULL;
3563 atom *OtherWalker = NULL;
3564 double *matrix = ReturnFullMatrixforSymmetric(cell_size);
3565 enum Shading *ColorList = NULL;
3566 double tmp;
3567 vector TranslationVector;
3568 //class StackClass<atom *> *CompStack = NULL;
3569 class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
3570 bool flag = true;
3571
3572// *out << Verbose(1) << "Begin of ScanForPeriodicCorrection." << endl;
3573
3574 ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");
3575 while (flag) {
3576 // remove bonds that are beyond bonddistance
3577 for(int i=0;i<NDIM;i++)
3578 TranslationVector.x[i] = 0.;
3579 // scan all bonds
3580 Binder = first;
3581 flag = false;
3582 while ((!flag) && (Binder->next != last)) {
3583 Binder = Binder->next;
3584 for (int i=0;i<NDIM;i++) {
3585 tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
3586 //*out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
3587 if (tmp > BondDistance) {
3588 OtherBinder = Binder->next; // note down binding partner for later re-insertion
3589 unlink(Binder); // unlink bond
3590// *out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
3591 flag = true;
3592 break;
3593 }
3594 }
3595 }
3596 if (flag) {
3597 // create translation vector from their periodically modified distance
3598 for (int i=0;i<NDIM;i++) {
3599 tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
3600 if (fabs(tmp) > BondDistance)
3601 TranslationVector.x[i] = (tmp < 0) ? +1. : -1.;
3602 }
3603 TranslationVector.MatrixMultiplication(matrix);
3604 //*out << "Translation vector is ";
3605 //TranslationVector.Output(out);
3606 //*out << endl;
3607 // apply to all atoms of first component via BFS
3608 for (int i=0;i<AtomCount;i++)
3609 ColorList[i] = white;
3610 AtomStack->Push(Binder->leftatom);
3611 while (!AtomStack->IsEmpty()) {
3612 Walker = AtomStack->PopFirst();
3613 //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
3614 ColorList[Walker->nr] = black; // mark as explored
3615 Walker->x.AddVector(&TranslationVector); // translate
3616 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through all binding partners
3617 if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
3618 OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
3619 if (ColorList[OtherWalker->nr] == white) {
3620 AtomStack->Push(OtherWalker); // push if yet unexplored
3621 }
3622 }
3623 }
3624 }
3625 // re-add bond
3626 link(Binder, OtherBinder);
3627 } else {
3628// *out << Verbose(2) << "No corrections for this fragment." << endl;
3629 }
3630 //delete(CompStack);
3631 }
3632
3633 // free allocated space from ReturnFullMatrixforSymmetric()
3634 delete(AtomStack);
3635 Free((void **)&ColorList, "molecule::ScanForPeriodicCorrection: *ColorList");
3636 Free((void **)&matrix, "molecule::ScanForPeriodicCorrection: *matrix");
3637// *out << Verbose(1) << "End of ScanForPeriodicCorrection." << endl;
3638};
3639
3640/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
3641 * \param *symm 6-dim array of unique symmetric matrix components
3642 * \return allocated NDIM*NDIM array with the symmetric matrix
3643 */
3644double * molecule::ReturnFullMatrixforSymmetric(double *symm)
3645{
3646 double *matrix = (double *) Malloc(sizeof(double)*NDIM*NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
3647 matrix[0] = symm[0];
3648 matrix[1] = symm[1];
3649 matrix[2] = symm[3];
3650 matrix[3] = symm[1];
3651 matrix[4] = symm[2];
3652 matrix[5] = symm[4];
3653 matrix[6] = symm[3];
3654 matrix[7] = symm[4];
3655 matrix[8] = symm[5];
3656 return matrix;
3657};
3658
3659bool KeyCompare::operator() (const KeySet SubgraphA, const KeySet SubgraphB) const
3660{
3661 //cout << "my check is used." << endl;
3662 if (SubgraphA.size() < SubgraphB.size()) {
3663 return true;
3664 } else {
3665 if (SubgraphA.size() > SubgraphB.size()) {
3666 return false;
3667 } else {
3668 KeySet::iterator IteratorA = SubgraphA.begin();
3669 KeySet::iterator IteratorB = SubgraphB.begin();
3670 while ((IteratorA != SubgraphA.end()) && (IteratorB != SubgraphB.end())) {
3671 if ((*IteratorA) < (*IteratorB))
3672 return true;
3673 else if ((*IteratorA) > (*IteratorB)) {
3674 return false;
3675 } // else, go on to next index
3676 IteratorA++;
3677 IteratorB++;
3678 } // end of while loop
3679 }// end of check in case of equal sizes
3680 }
3681 return false; // if we reach this point, they are equal
3682};
3683
3684//bool operator < (KeySet SubgraphA, KeySet SubgraphB)
3685//{
3686// return KeyCompare(SubgraphA, SubgraphB);
3687//};
3688
3689/** Checking whether KeySet is not already present in Graph, if so just adds factor.
3690 * \param *out output stream for debugging
3691 * \param &set KeySet to insert
3692 * \param &graph Graph to insert into
3693 * \param *counter pointer to unique fragment count
3694 * \param factor energy factor for the fragment
3695 */
3696inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment)
3697{
3698 GraphTestPair testGraphInsert;
3699
3700 testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,Fragment->TEFactor))); // store fragment number and current factor
3701 if (testGraphInsert.second) {
3702 *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " successfully inserted." << endl;
3703 Fragment->FragmentCounter++;
3704 } else {
3705 *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
3706 ((*(testGraphInsert.first)).second).second += Fragment->TEFactor; // increase the "created" counter
3707 *out << Verbose(2) << "New factor is " << ((*(testGraphInsert.first)).second).second << "." << endl;
3708 }
3709};
3710//void inline InsertIntoGraph(ofstream *out, KeyStack &stack, Graph &graph, int *counter, double factor)
3711//{
3712// // copy stack contents to set and call overloaded function again
3713// KeySet set;
3714// for(KeyStack::iterator runner = stack.begin(); runner != stack.begin(); runner++)
3715// set.insert((*runner));
3716// InsertIntoGraph(out, set, graph, counter, factor);
3717//};
3718
3719/** Inserts each KeySet in \a graph2 into \a graph1.
3720 * \param *out output stream for debugging
3721 * \param graph1 first (dest) graph
3722 * \param graph2 second (source) graph
3723 * \param *counter keyset counter that gets increased
3724 */
3725inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
3726{
3727 GraphTestPair testGraphInsert;
3728
3729 for(Graph::iterator runner = graph2.begin(); runner != graph2.end(); runner++) {
3730 testGraphInsert = graph1.insert(GraphPair ((*runner).first,pair<int,double>((*counter)++,((*runner).second).second))); // store fragment number and current factor
3731 if (testGraphInsert.second) {
3732 *out << Verbose(2) << "KeySet " << (*counter)-1 << " successfully inserted." << endl;
3733 } else {
3734 *out << Verbose(2) << "KeySet " << (*counter)-1 << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
3735 ((*(testGraphInsert.first)).second).second += (*runner).second.second;
3736 *out << Verbose(2) << "New factor is " << (*(testGraphInsert.first)).second.second << "." << endl;
3737 }
3738 }
3739};
3740
3741
3742/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
3743 * -# constructs a complete keyset of the molecule
3744 * -# In a loop over all possible roots from the given rootstack
3745 * -# increases order of root site
3746 * -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
3747 * -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
3748as the restricted one and each site in the set as the root)
3749 * -# these are merged into a fragment list of keysets
3750 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
3751 * Important only is that we create all fragments, it is not important if we create them more than once
3752 * as these copies are filtered out via use of the hash table (KeySet).
3753 * \param *out output stream for debugging
3754 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
3755 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
3756 * \param *MinimumRingSize minimum ring size for each atom (molecule::Atomcount)
3757 * \return pointer to Graph list
3758 */
3759void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
3760{
3761 Graph ***FragmentLowerOrdersList = NULL;
3762 int Order, NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
3763 int counter = 0;
3764 int UpgradeCount = RootStack.size();
3765 KeyStack FragmentRootStack;
3766 int RootKeyNr, RootNr;
3767 struct UniqueFragments FragmentSearch;
3768
3769 *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
3770
3771 // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
3772 // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
3773 NumMoleculesOfOrder = (int *) Malloc(sizeof(int)*UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
3774 FragmentLowerOrdersList = (Graph ***) Malloc(sizeof(Graph **)*UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
3775
3776 // initialise the fragments structure
3777 FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
3778 FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
3779 FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
3780 FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
3781 FragmentSearch.FragmentCounter = 0;
3782 FragmentSearch.FragmentSet = new KeySet;
3783 FragmentSearch.Root = FindAtom(RootKeyNr);
3784 for (int i=0;i<AtomCount;i++) {
3785 FragmentSearch.Labels[i] = -1;
3786 FragmentSearch.ShortestPathList[i] = -1;
3787 }
3788
3789 // Construct the complete KeySet which we need for topmost level only (but for all Roots)
3790 atom *Walker = start;
3791 KeySet CompleteMolecule;
3792 while (Walker->next != end) {
3793 Walker = Walker->next;
3794 CompleteMolecule.insert(Walker->GetTrueFather()->nr);
3795 }
3796
3797 // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
3798 // each has to be split up. E.g. for the second level we have one from 5th, one from 4th, two from 3th (which in turn is one from 5th, one from 4th),
3799 // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
3800 // with the order along the cells as this: 5433222211111111 for BondOrder 5 needing 16=pow(2,5-1) cells (only we use bit-shifting which is faster)
3801 RootNr = 0; // counts through the roots in RootStack
3802 while (RootNr < UpgradeCount) {
3803 RootKeyNr = RootStack.front();
3804 RootStack.pop_front();
3805 Walker = FindAtom(RootKeyNr);
3806 // check cyclic lengths
3807 if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 >= MinimumRingSize[Walker->GetTrueFather()->nr])) {
3808 *out << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
3809 } else {
3810 // increase adaptive order by one
3811 Walker->GetTrueFather()->AdaptiveOrder++;
3812 Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
3813
3814 // initialise Order-dependent entries of UniqueFragments structure
3815 FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
3816 FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
3817 for (int i=0;i<Order;i++) {
3818 FragmentSearch.BondsPerSPList[2*i] = new bond(); // start node
3819 FragmentSearch.BondsPerSPList[2*i+1] = new bond(); // end node
3820 FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1]; // intertwine these two
3821 FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
3822 FragmentSearch.BondsPerSPCount[i] = 0;
3823 }
3824
3825 // allocate memory for all lower level orders in this 1D-array of ptrs
3826 NumLevels = 1 << (Order-1); // (int)pow(2,Order);
3827 FragmentLowerOrdersList[RootNr] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
3828
3829 // create top order where nothing is reduced
3830 *out << Verbose(0) << "==============================================================================================================" << endl;
3831 *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr-1) << " Roots remaining." << endl;
3832
3833 // Create list of Graphs of current Bond Order (i.e. F_{ij})
3834 FragmentLowerOrdersList[RootNr][0] = new Graph;
3835 FragmentSearch.TEFactor = 1.;
3836 FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0]; // set to insertion graph
3837 FragmentSearch.Root = Walker;
3838 NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
3839 *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
3840 NumMolecules = 0;
3841
3842 if ((NumLevels >> 1) > 0) {
3843 // create lower order fragments
3844 *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
3845 Order = Walker->AdaptiveOrder;
3846 for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
3847 // step down to next order at (virtual) boundary of powers of 2 in array
3848 while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
3849 Order--;
3850 *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
3851 for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
3852 int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
3853 *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
3854 *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
3855
3856 // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
3857 //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
3858 //NumMolecules = 0;
3859 FragmentLowerOrdersList[RootNr][dest] = new Graph;
3860 for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
3861 for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
3862 Graph TempFragmentList;
3863 FragmentSearch.TEFactor = -(*runner).second.second;
3864 FragmentSearch.Leaflet = &TempFragmentList; // set to insertion graph
3865 FragmentSearch.Root = FindAtom(*sprinter);
3866 NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
3867 // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
3868 *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
3869 InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
3870 }
3871 }
3872 *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
3873 }
3874 }
3875 }
3876 // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
3877 //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
3878 TotalNumMolecules += NumMoleculesOfOrder[RootNr];
3879 *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
3880 RootStack.push_back(RootKeyNr); // put back on stack
3881 RootNr++;
3882
3883 // free Order-dependent entries of UniqueFragments structure for next loop cycle
3884 Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
3885 for (int i=0;i<Order;i++) {
3886 delete(FragmentSearch.BondsPerSPList[2*i]);
3887 delete(FragmentSearch.BondsPerSPList[2*i+1]);
3888 }
3889 Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
3890 }
3891 }
3892 *out << Verbose(0) << "==============================================================================================================" << endl;
3893 *out << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl;
3894 *out << Verbose(0) << "==============================================================================================================" << endl;
3895
3896 // cleanup FragmentSearch structure
3897 Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
3898 Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
3899 delete(FragmentSearch.FragmentSet);
3900
3901 // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
3902 // 5433222211111111
3903 // 43221111
3904 // 3211
3905 // 21
3906 // 1
3907
3908 // Subsequently, we combine all into a single list (FragmentList)
3909
3910 *out << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl;
3911 if (FragmentList == NULL) {
3912 FragmentList = new Graph;
3913 counter = 0;
3914 } else {
3915 counter = FragmentList->size();
3916 }
3917 RootNr = 0;
3918 while (!RootStack.empty()) {
3919 RootKeyNr = RootStack.front();
3920 RootStack.pop_front();
3921 Walker = FindAtom(RootKeyNr);
3922 NumLevels = 1 << (Walker->AdaptiveOrder - 1);
3923 for(int i=0;i<NumLevels;i++) {
3924 InsertGraphIntoGraph(out, *FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
3925 delete(FragmentLowerOrdersList[RootNr][i]);
3926 }
3927 Free((void **)&FragmentLowerOrdersList[RootNr], "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
3928 RootNr++;
3929 }
3930 Free((void **)&FragmentLowerOrdersList, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
3931 Free((void **)&NumMoleculesOfOrder, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
3932
3933 *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
3934};
3935
3936/** Comparison function for GSL heapsort on distances in two molecules.
3937 * \param *a
3938 * \param *b
3939 * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
3940 */
3941int CompareDoubles (const void * a, const void * b)
3942{
3943 if (*(double *)a > *(double *)b)
3944 return -1;
3945 else if (*(double *)a < *(double *)b)
3946 return 1;
3947 else
3948 return 0;
3949};
3950
3951/** Determines whether two molecules actually contain the same atoms and coordination.
3952 * \param *out output stream for debugging
3953 * \param *OtherMolecule the molecule to compare this one to
3954 * \param threshold upper limit of difference when comparing the coordination.
3955 * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
3956 */
3957int * molecule::IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold)
3958{
3959 int flag;
3960 double *Distances = NULL, *OtherDistances = NULL;
3961 vector CenterOfGravity, OtherCenterOfGravity;
3962 size_t *PermMap = NULL, *OtherPermMap = NULL;
3963 int *PermutationMap = NULL;
3964 atom *Walker = NULL;
3965 bool result = true; // status of comparison
3966
3967 *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
3968 /// first count both their atoms and elements and update lists thereby ...
3969 //*out << Verbose(0) << "Counting atoms, updating list" << endl;
3970 CountAtoms(out);
3971 OtherMolecule->CountAtoms(out);
3972 CountElements();
3973 OtherMolecule->CountElements();
3974
3975 /// ... and compare:
3976 /// -# AtomCount
3977 if (result) {
3978 if (AtomCount != OtherMolecule->AtomCount) {
3979 *out << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
3980 result = false;
3981 } else *out << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
3982 }
3983 /// -# ElementCount
3984 if (result) {
3985 if (ElementCount != OtherMolecule->ElementCount) {
3986 *out << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
3987 result = false;
3988 } else *out << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
3989 }
3990 /// -# ElementsInMolecule
3991 if (result) {
3992 for (flag=0;flag<MAX_ELEMENTS;flag++) {
3993 //*out << Verbose(5) << "Element " << flag << ": " << ElementsInMolecule[flag] << " <-> " << OtherMolecule->ElementsInMolecule[flag] << "." << endl;
3994 if (ElementsInMolecule[flag] != OtherMolecule->ElementsInMolecule[flag])
3995 break;
3996 }
3997 if (flag < MAX_ELEMENTS) {
3998 *out << Verbose(4) << "ElementsInMolecule don't match." << endl;
3999 result = false;
4000 } else *out << Verbose(4) << "ElementsInMolecule match." << endl;
4001 }
4002 /// then determine and compare center of gravity for each molecule ...
4003 if (result) {
4004 *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
4005 DetermineCenterOfGravity(CenterOfGravity);
4006 OtherMolecule->DetermineCenterOfGravity(OtherCenterOfGravity);
4007 *out << Verbose(5) << "Center of Gravity: ";
4008 CenterOfGravity.Output(out);
4009 *out << endl << Verbose(5) << "Other Center of Gravity: ";
4010 OtherCenterOfGravity.Output(out);
4011 *out << endl;
4012 if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
4013 *out << Verbose(4) << "Centers of gravity don't match." << endl;
4014 result = false;
4015 }
4016 }
4017
4018 /// ... then make a list with the euclidian distance to this center for each atom of both molecules
4019 if (result) {
4020 *out << Verbose(5) << "Calculating distances" << endl;
4021 Distances = (double *) Malloc(sizeof(double)*AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
4022 OtherDistances = (double *) Malloc(sizeof(double)*AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
4023 Walker = start;
4024 while (Walker->next != end) {
4025 Walker = Walker->next;
4026 //for (i=0;i<AtomCount;i++) {
4027 Distances[Walker->nr] = CenterOfGravity.Distance(&Walker->x);
4028 }
4029 Walker = OtherMolecule->start;
4030 while (Walker->next != OtherMolecule->end) {
4031 Walker = Walker->next;
4032 OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
4033 }
4034
4035 /// ... sort each list (using heapsort (o(N log N)) from GSL)
4036 *out << Verbose(5) << "Sorting distances" << endl;
4037 PermMap = (size_t *) Malloc(sizeof(size_t)*AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
4038 OtherPermMap = (size_t *) Malloc(sizeof(size_t)*AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
4039 gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
4040 gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
4041 PermutationMap = (int *) Malloc(sizeof(int)*AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
4042 *out << Verbose(5) << "Combining Permutation Maps" << endl;
4043 for(int i=0;i<AtomCount;i++)
4044 PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
4045
4046 /// ... and compare them step by step, whether the difference is individiually(!) below \a threshold for all
4047 *out << Verbose(4) << "Comparing distances" << endl;
4048 flag = 0;
4049 for (int i=0;i<AtomCount;i++) {
4050 *out << Verbose(5) << "Distances: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl;
4051 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold)
4052 flag = 1;
4053 }
4054 Free((void **)&PermMap, "molecule::IsEqualToWithinThreshold: *PermMap");
4055 Free((void **)&OtherPermMap, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
4056
4057 /// free memory
4058 Free((void **)&Distances, "molecule::IsEqualToWithinThreshold: Distances");
4059 Free((void **)&OtherDistances, "molecule::IsEqualToWithinThreshold: OtherDistances");
4060 if (flag) { // if not equal
4061 Free((void **)&PermutationMap, "molecule::IsEqualToWithinThreshold: *PermutationMap");
4062 result = false;
4063 }
4064 }
4065 /// return pointer to map if all distances were below \a threshold
4066 *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
4067 if (result) {
4068 *out << Verbose(3) << "Result: Equal." << endl;
4069 return PermutationMap;
4070 } else {
4071 *out << Verbose(3) << "Result: Not equal." << endl;
4072 return NULL;
4073 }
4074};
4075
4076/** Returns an index map for two father-son-molecules.
4077 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
4078 * \param *out output stream for debugging
4079 * \param *OtherMolecule corresponding molecule with fathers
4080 * \return allocated map of size molecule::AtomCount with map
4081 * \todo make this with a good sort O(n), not O(n^2)
4082 */
4083int * molecule::GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule)
4084{
4085 atom *Walker = NULL, *OtherWalker = NULL;
4086 *out << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
4087 int *AtomicMap = (int *) Malloc(sizeof(int)*AtomCount, "molecule::GetAtomicMap: *AtomicMap"); //Calloc
4088 for (int i=0;i<AtomCount;i++)
4089 AtomicMap[i] = -1;
4090 if (OtherMolecule == this) { // same molecule
4091 for (int i=0;i<AtomCount;i++) // no need as -1 means already that there is trivial correspondence
4092 AtomicMap[i] = i;
4093 *out << Verbose(4) << "Map is trivial." << endl;
4094 } else {
4095 *out << Verbose(4) << "Map is ";
4096 Walker = start;
4097 while (Walker->next != end) {
4098 Walker = Walker->next;
4099 if (Walker->father == NULL) {
4100 AtomicMap[Walker->nr] = -2;
4101 } else {
4102 OtherWalker = OtherMolecule->start;
4103 while (OtherWalker->next != OtherMolecule->end) {
4104 OtherWalker = OtherWalker->next;
4105 //for (int i=0;i<AtomCount;i++) { // search atom
4106 //for (int j=0;j<OtherMolecule->AtomCount;j++) {
4107 //*out << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
4108 if (Walker->father == OtherWalker)
4109 AtomicMap[Walker->nr] = OtherWalker->nr;
4110 }
4111 }
4112 *out << AtomicMap[Walker->nr] << "\t";
4113 }
4114 *out << endl;
4115 }
4116 *out << Verbose(3) << "End of GetFatherAtomicMap." << endl;
4117 return AtomicMap;
4118};
4119
Note: See TracBrowser for help on using the repository browser.