source: molecuilder/src/molecules.cpp@ 747b10

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

CheckOrderAtSite: BUGFIX - Added MinimumRingSize to parameter list, as otherwise infinite loop

If FragmentBOSSANOVA checks for remaining below MinimumRingSize per atom, CheckOrderAtSite however does not and fills AtomMask with values for atoms not to be included, we end up with an endless loop. Hence, MinimumRingSize is given as an additional parameter and CheckOrderAtSite checks for it for adaptive and for global increase (for stepwise we don't have the loop problem and FragmentBOSSANOVA will do the right thing).

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