source: molecuilder/src/molecules.cpp@ 5eb05a

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

molecule::CreateAdjacencyList() now needs IsAngstroem as parameter

This is necessary, as the database values (covalent radii et al) are in Angstroem, hence don't match if we use Bohrradii as length unit instead. CreateAdjacencyList() converts units and now finds correct bond structure and fragments properly.

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