source: molecuilder/src/molecules.cpp@ 3f805d

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

CyclicStructureAnalysis: BUGFIX - MinimumRingSize was not set correctly

  • During BFS search from non-loop members to loops, Walker was pushed instead of Root at start of BFS
  • if Backedge lead to a terminal hydrogen, this hydrogen ended up with AtomCount MinimumRingSize instead of correct one.
  • Property mode set to 100644
File size: 191.1 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 MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
1891 while (Walker != BackEdge->rightatom) {
1892 Walker = PredecessorList[Walker->nr];
1893 if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
1894 MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
1895 }
1896 if ((RingSize < MinRingSize) || (MinRingSize == -1))
1897 MinRingSize = RingSize;
1898 } else {
1899 *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
1900 }
1901
1902 // now clean the lists
1903 while (!TouchedStack->IsEmpty()){
1904 Walker = TouchedStack->PopFirst();
1905 PredecessorList[Walker->nr] = NULL;
1906 ShortestPathList[Walker->nr] = -1;
1907 ColorList[Walker->nr] = white;
1908 }
1909 }
1910 if (MinRingSize != -1) {
1911 // go over all atoms
1912 Root = start;
1913 while(Root->next != end) {
1914 Root = Root->next;
1915
1916 if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
1917 Walker = Root;
1918 ShortestPathList[Walker->nr] = 0;
1919 BFSStack->ClearStack(); // start with empty BFS stack
1920 BFSStack->Push(Walker);
1921 TouchedStack->Push(Walker);
1922 //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
1923 OtherAtom = Walker;
1924 while (OtherAtom != NULL) { // look for Root
1925 Walker = BFSStack->PopFirst();
1926 //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
1927 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
1928 Binder = ListOfBondsPerAtom[Walker->nr][i];
1929 if ((Binder != BackEdge) || (NumberOfBondsPerAtom[Walker->nr] == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
1930 OtherAtom = Binder->GetOtherAtom(Walker);
1931 //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
1932 if (ColorList[OtherAtom->nr] == white) {
1933 TouchedStack->Push(OtherAtom);
1934 ColorList[OtherAtom->nr] = lightgray;
1935 PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
1936 ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
1937 //*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;
1938 if (MinimumRingSize[OtherAtom->GetTrueFather()->nr] != AtomCount) { // if the other atom is connected to a ring
1939 MinimumRingSize[Root->GetTrueFather()->nr] = ShortestPathList[OtherAtom->nr]+MinimumRingSize[OtherAtom->GetTrueFather()->nr];
1940 OtherAtom = NULL; //break;
1941 break;
1942 } else
1943 BFSStack->Push(OtherAtom);
1944 } else {
1945 //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
1946 }
1947 } else {
1948 //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
1949 }
1950 }
1951 ColorList[Walker->nr] = black;
1952 //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
1953 }
1954
1955 // now clean the lists
1956 while (!TouchedStack->IsEmpty()){
1957 Walker = TouchedStack->PopFirst();
1958 PredecessorList[Walker->nr] = NULL;
1959 ShortestPathList[Walker->nr] = -1;
1960 ColorList[Walker->nr] = white;
1961 }
1962 }
1963 *out << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl;
1964 }
1965 *out << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl;
1966 } else
1967 *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
1968
1969 Free((void **)&PredecessorList, "molecule::CyclicStructureAnalysis: **PredecessorList");
1970 Free((void **)&ShortestPathList, "molecule::CyclicStructureAnalysis: **ShortestPathList");
1971 Free((void **)&ColorList, "molecule::CyclicStructureAnalysis: **ColorList");
1972 delete(BFSStack);
1973};
1974
1975/** Sets the next component number.
1976 * This is O(N) as the number of bonds per atom is bound.
1977 * \param *vertex atom whose next atom::*ComponentNr is to be set
1978 * \param nr number to use
1979 */
1980void molecule::SetNextComponentNumber(atom *vertex, int nr)
1981{
1982 int i=0;
1983 if (vertex != NULL) {
1984 for(;i<NumberOfBondsPerAtom[vertex->nr];i++) {
1985 if (vertex->ComponentNr[i] == -1) { // check if not yet used
1986 vertex->ComponentNr[i] = nr;
1987 break;
1988 }
1989 else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
1990 break; // breaking here will not cause error!
1991 }
1992 if (i == NumberOfBondsPerAtom[vertex->nr])
1993 cerr << "Error: All Component entries are already occupied!" << endl;
1994 } else
1995 cerr << "Error: Given vertex is NULL!" << endl;
1996};
1997
1998/** Output a list of flags, stating whether the bond was visited or not.
1999 * \param *out output stream for debugging
2000 */
2001void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
2002{
2003 for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
2004 *out << vertex->ComponentNr[i] << " ";
2005};
2006
2007/** Allocates memory for all atom::*ComponentNr in this molecule and sets each entry to -1.
2008 */
2009void molecule::InitComponentNumbers()
2010{
2011 atom *Walker = start;
2012 while(Walker->next != end) {
2013 Walker = Walker->next;
2014 if (Walker->ComponentNr != NULL)
2015 Free((void **)&Walker->ComponentNr, "molecule::InitComponentNumbers: **Walker->ComponentNr");
2016 Walker->ComponentNr = (int *) Malloc(sizeof(int)*NumberOfBondsPerAtom[Walker->nr], "molecule::InitComponentNumbers: *Walker->ComponentNr");
2017 for (int i=NumberOfBondsPerAtom[Walker->nr];i--;)
2018 Walker->ComponentNr[i] = -1;
2019 }
2020};
2021
2022/** Returns next unused bond for this atom \a *vertex or NULL of none exists.
2023 * \param *vertex atom to regard
2024 * \return bond class or NULL
2025 */
2026bond * molecule::FindNextUnused(atom *vertex)
2027{
2028 for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
2029 if (ListOfBondsPerAtom[vertex->nr][i]->IsUsed() == white)
2030 return(ListOfBondsPerAtom[vertex->nr][i]);
2031 return NULL;
2032};
2033
2034/** Resets bond::Used flag of all bonds in this molecule.
2035 * \return true - success, false - -failure
2036 */
2037void molecule::ResetAllBondsToUnused()
2038{
2039 bond *Binder = first;
2040 while (Binder->next != last) {
2041 Binder = Binder->next;
2042 Binder->ResetUsed();
2043 }
2044};
2045
2046/** Resets atom::nr to -1 of all atoms in this molecule.
2047 */
2048void molecule::ResetAllAtomNumbers()
2049{
2050 atom *Walker = start;
2051 while (Walker->next != end) {
2052 Walker = Walker->next;
2053 Walker->GraphNr = -1;
2054 }
2055};
2056
2057/** Output a list of flags, stating whether the bond was visited or not.
2058 * \param *out output stream for debugging
2059 * \param *list
2060 */
2061void OutputAlreadyVisited(ofstream *out, int *list)
2062{
2063 *out << Verbose(4) << "Already Visited Bonds:\t";
2064 for(int i=1;i<=list[0];i++) *out << Verbose(0) << list[i] << " ";
2065 *out << endl;
2066};
2067
2068/** Estimates by educated guessing (using upper limit) the expected number of fragments.
2069 * The upper limit is
2070 * \f[
2071 * n = N \cdot C^k
2072 * \f]
2073 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
2074 * \param *out output stream for debugging
2075 * \param order bond order k
2076 * \return number n of fragments
2077 */
2078int molecule::GuesstimateFragmentCount(ofstream *out, int order)
2079{
2080 int c = 0;
2081 int FragmentCount;
2082 // get maximum bond degree
2083 atom *Walker = start;
2084 while (Walker->next != end) {
2085 Walker = Walker->next;
2086 c = (NumberOfBondsPerAtom[Walker->nr] > c) ? NumberOfBondsPerAtom[Walker->nr] : c;
2087 }
2088 FragmentCount = NoNonHydrogen*(1 << (c*order));
2089 *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
2090 return FragmentCount;
2091};
2092
2093/** Scans a single line for number and puts them into \a KeySet.
2094 * \param *out output stream for debugging
2095 * \param *buffer buffer to scan
2096 * \param &CurrentSet filled KeySet on return
2097 * \return true - at least one valid atom id parsed, false - CurrentSet is empty
2098 */
2099bool molecule::ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
2100{
2101 stringstream line;
2102 int AtomNr;
2103 int status = 0;
2104
2105 line.str(buffer);
2106 while (!line.eof()) {
2107 line >> AtomNr;
2108 if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
2109 CurrentSet.insert(AtomNr); // insert at end, hence in same order as in file!
2110 status++;
2111 } // else it's "-1" or else and thus must not be added
2112 }
2113 *out << Verbose(1) << "The scanned KeySet is ";
2114 for(KeySet::iterator runner = CurrentSet.begin(); runner != CurrentSet.end(); runner++) {
2115 *out << (*runner) << "\t";
2116 }
2117 *out << endl;
2118 return (status != 0);
2119};
2120
2121/** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
2122 * Does two-pass scanning:
2123 * -# Scans the keyset file and initialises a temporary graph
2124 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
2125 * Finally, the temporary graph is inserted into the given \a FragmentList for return.
2126 * \param *out output stream for debugging
2127 * \param *path path to file
2128 * \param *FragmentList empty, filled on return
2129 * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
2130 */
2131bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
2132{
2133 bool status = true;
2134 ifstream InputFile;
2135 stringstream line;
2136 GraphTestPair testGraphInsert;
2137 int NumberOfFragments = 0;
2138 double TEFactor;
2139 char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
2140
2141 if (FragmentList == NULL) { // check list pointer
2142 FragmentList = new Graph;
2143 }
2144
2145 // 1st pass: open file and read
2146 *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
2147 sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
2148 InputFile.open(filename);
2149 if (InputFile != NULL) {
2150 // each line represents a new fragment
2151 char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - *buffer");
2152 // 1. parse keysets and insert into temp. graph
2153 while (!InputFile.eof()) {
2154 InputFile.getline(buffer, MAXSTRINGSIZE);
2155 KeySet CurrentSet;
2156 if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(out, buffer, CurrentSet))) { // if at least one valid atom was added, write config
2157 testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1))); // store fragment number and current factor
2158 if (!testGraphInsert.second) {
2159 cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
2160 }
2161 //FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
2162 }
2163 }
2164 // 2. Free and done
2165 InputFile.close();
2166 InputFile.clear();
2167 Free((void **)&buffer, "molecule::ParseKeySetFile - *buffer");
2168 *out << Verbose(1) << "done." << endl;
2169 } else {
2170 *out << Verbose(1) << "File " << filename << " not found." << endl;
2171 status = false;
2172 }
2173
2174 // 2nd pass: open TEFactors file and read
2175 *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
2176 sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, TEFACTORSFILE);
2177 InputFile.open(filename);
2178 if (InputFile != NULL) {
2179 // 3. add found TEFactors to each keyset
2180 NumberOfFragments = 0;
2181 for(Graph::iterator runner = FragmentList->begin();runner != FragmentList->end(); runner++) {
2182 if (!InputFile.eof()) {
2183 InputFile >> TEFactor;
2184 (*runner).second.second = TEFactor;
2185 *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
2186 } else {
2187 status = false;
2188 break;
2189 }
2190 }
2191 // 4. Free and done
2192 InputFile.close();
2193 *out << Verbose(1) << "done." << endl;
2194 } else {
2195 *out << Verbose(1) << "File " << filename << " not found." << endl;
2196 status = false;
2197 }
2198
2199 // free memory
2200 Free((void **)&filename, "molecule::ParseKeySetFile - filename");
2201
2202 return status;
2203};
2204
2205/** Stores keysets and TEFactors to file.
2206 * \param *out output stream for debugging
2207 * \param KeySetList Graph with Keysets and factors
2208 * \param *path path to file
2209 * \return true - file written successfully, false - writing failed
2210 */
2211bool molecule::StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
2212{
2213 ofstream output;
2214 bool status = true;
2215 string line;
2216
2217 // open KeySet file
2218 line = path;
2219 line.append("/");
2220 line += FRAGMENTPREFIX;
2221 line += KEYSETFILE;
2222 output.open(line.c_str(), ios::out);
2223 *out << Verbose(1) << "Saving key sets of the total graph ... ";
2224 if(output != NULL) {
2225 for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
2226 for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
2227 if (sprinter != (*runner).first.begin())
2228 output << "\t";
2229 output << *sprinter;
2230 }
2231 output << endl;
2232 }
2233 *out << "done." << endl;
2234 } else {
2235 cerr << "Unable to open " << line << " for writing keysets!" << endl;
2236 status = false;
2237 }
2238 output.close();
2239 output.clear();
2240
2241 // open TEFactors file
2242 line = path;
2243 line.append("/");
2244 line += FRAGMENTPREFIX;
2245 line += TEFACTORSFILE;
2246 output.open(line.c_str(), ios::out);
2247 *out << Verbose(1) << "Saving TEFactors of the total graph ... ";
2248 if(output != NULL) {
2249 for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++)
2250 output << (*runner).second.second << endl;
2251 *out << Verbose(1) << "done." << endl;
2252 } else {
2253 *out << Verbose(1) << "failed to open " << line << "." << endl;
2254 status = false;
2255 }
2256 output.close();
2257
2258 return status;
2259};
2260
2261/** Storing the bond structure of a molecule to file.
2262 * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
2263 * \param *out output stream for debugging
2264 * \param *path path to file
2265 * \return true - file written successfully, false - writing failed
2266 */
2267bool molecule::StoreAdjacencyToFile(ofstream *out, char *path)
2268{
2269 ofstream AdjacencyFile;
2270 atom *Walker = NULL;
2271 stringstream line;
2272 bool status = true;
2273
2274 line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
2275 AdjacencyFile.open(line.str().c_str(), ios::out);
2276 *out << Verbose(1) << "Saving adjacency list ... ";
2277 if (AdjacencyFile != NULL) {
2278 Walker = start;
2279 while(Walker->next != end) {
2280 Walker = Walker->next;
2281 AdjacencyFile << Walker->nr << "\t";
2282 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
2283 AdjacencyFile << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
2284 AdjacencyFile << endl;
2285 }
2286 AdjacencyFile.close();
2287 *out << Verbose(1) << "done." << endl;
2288 } else {
2289 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
2290 status = false;
2291 }
2292
2293 return status;
2294};
2295
2296/** Checks contents of adjacency file against bond structure in structure molecule.
2297 * \param *out output stream for debugging
2298 * \param *path path to file
2299 * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
2300 * \return true - structure is equal, false - not equivalence
2301 */
2302bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
2303{
2304 ifstream File;
2305 stringstream filename;
2306 bool status = true;
2307 char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
2308
2309 filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
2310 File.open(filename.str().c_str(), ios::out);
2311 *out << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... ";
2312 if (File != NULL) {
2313 // allocate storage structure
2314 int NonMatchNumber = 0; // will number of atoms with differing bond structure
2315 int *CurrentBonds = (int *) Malloc(sizeof(int)*8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
2316 int CurrentBondsOfAtom;
2317
2318 // Parse the file line by line and count the bonds
2319 while (!File.eof()) {
2320 File.getline(buffer, MAXSTRINGSIZE);
2321 stringstream line;
2322 line.str(buffer);
2323 int AtomNr = -1;
2324 line >> AtomNr;
2325 CurrentBondsOfAtom = -1; // we count one too far due to line end
2326 // parse into structure
2327 if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
2328 while (!line.eof())
2329 line >> CurrentBonds[ ++CurrentBondsOfAtom ];
2330 // compare against present bonds
2331 //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
2332 if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
2333 for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
2334 int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
2335 int j = 0;
2336 for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
2337 if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
2338 ListOfAtoms[AtomNr] = NULL;
2339 NonMatchNumber++;
2340 status = false;
2341 //out << "[" << id << "]\t";
2342 } else {
2343 //out << id << "\t";
2344 }
2345 }
2346 //out << endl;
2347 } else {
2348 *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
2349 status = false;
2350 }
2351 }
2352 }
2353 File.close();
2354 File.clear();
2355 if (status) { // if equal we parse the KeySetFile
2356 *out << Verbose(1) << "done: Equal." << endl;
2357 status = true;
2358 } else
2359 *out << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl;
2360 Free((void **)&CurrentBonds, "molecule::CheckAdjacencyFileAgainstMolecule - **CurrentBonds");
2361 } else {
2362 *out << Verbose(1) << "Adjacency file not found." << endl;
2363 status = false;
2364 }
2365 *out << endl;
2366 Free((void **)&buffer, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
2367
2368 return status;
2369};
2370
2371/** Checks whether the OrderAtSite is still below \a Order at some site.
2372 * \param *out output stream for debugging
2373 * \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
2374 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
2375 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
2376 * \param *MinimumRingSize array of max. possible order to avoid loops
2377 * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
2378 * \return true - needs further fragmentation, false - does not need fragmentation
2379 */
2380bool molecule::CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
2381{
2382 atom *Walker = start;
2383 bool status = false;
2384 ifstream InputFile;
2385
2386 // initialize mask list
2387 for(int i=AtomCount;i--;)
2388 AtomMask[i] = false;
2389
2390 if (Order < 0) { // adaptive increase of BondOrder per site
2391 if (AtomMask[AtomCount] == true) // break after one step
2392 return false;
2393 // parse the EnergyPerFragment file
2394 char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
2395 sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
2396 InputFile.open(buffer, ios::in);
2397 if ((InputFile != NULL) && (GlobalKeySetList != NULL)) {
2398 // transmorph graph keyset list into indexed KeySetList
2399 map<int,KeySet> IndexKeySetList;
2400 for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
2401 IndexKeySetList.insert( pair<int,KeySet>(runner->second.first,runner->first) );
2402 }
2403 int lines = 0;
2404 // count the number of lines, i.e. the number of fragments
2405 InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
2406 InputFile.getline(buffer, MAXSTRINGSIZE);
2407 while(!InputFile.eof()) {
2408 InputFile.getline(buffer, MAXSTRINGSIZE);
2409 lines++;
2410 }
2411 *out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl; // one endline too much
2412 InputFile.clear();
2413 InputFile.seekg(ios::beg);
2414 map<int, pair<double,int> > AdaptiveCriteriaList; // (Root No., (Value, Order)) !
2415 int No, FragOrder;
2416 double Value;
2417 // each line represents a fragment root (Atom::nr) id and its energy contribution
2418 InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
2419 InputFile.getline(buffer, MAXSTRINGSIZE);
2420 while(!InputFile.eof()) {
2421 InputFile.getline(buffer, MAXSTRINGSIZE);
2422 if (strlen(buffer) > 2) {
2423 //*out << Verbose(2) << "Scanning: " << buffer;
2424 stringstream line(buffer);
2425 line >> FragOrder;
2426 line >> ws >> No;
2427 line >> ws >> Value; // skip time entry
2428 line >> ws >> Value;
2429 No -= 1; // indices start at 1 in file, not 0
2430 //*out << Verbose(2) << " - yields (" << No << "," << Value << ")" << endl;
2431
2432 // clean the list of those entries that have been superceded by higher order terms already
2433 map<int,KeySet>::iterator marker = IndexKeySetList.find(No); // find keyset to Frag No.
2434 if (marker != IndexKeySetList.end()) { // if found
2435 // 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
2436 pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( Value, Order) ));
2437 map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
2438 if (!InsertedElement.second) { // this root is already present
2439 if ((*PresentItem).second.second < FragOrder) // if order there is lower, update entry with higher-order term
2440 //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)
2441 { // if value is smaller, update value and order
2442 (*PresentItem).second.first = Value;
2443 (*PresentItem).second.second = FragOrder;
2444 *out << Verbose(2) << "Updated element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
2445 }
2446 } else {
2447 *out << Verbose(2) << "Inserted element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
2448 }
2449 } else {
2450 *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
2451 }
2452 }
2453 }
2454 // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
2455 map<double, pair<int,int> > FinalRootCandidates;
2456 *out << Verbose(1) << "Root candidate list is: " << endl;
2457 for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList.begin(); runner != AdaptiveCriteriaList.end(); runner++) {
2458 Walker = FindAtom((*runner).first);
2459 if (Walker != NULL) {
2460 if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
2461 *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
2462 FinalRootCandidates.insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
2463 }
2464 } else {
2465 cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
2466 }
2467 }
2468 // pick the ones still below threshold and mark as to be adaptively updated
2469 for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
2470 No = (*runner).second.first;
2471 Walker = FindAtom(No);
2472 if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
2473 *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
2474 AtomMask[No] = true;
2475 status = true;
2476 } else
2477 *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;
2478 }
2479 // close and done
2480 InputFile.close();
2481 InputFile.clear();
2482 } else {
2483 cerr << "Unable to parse " << buffer << " file, incrementing all." << endl;
2484 while (Walker->next != end) {
2485 Walker = Walker->next;
2486 #ifdef ADDHYDROGEN
2487 if (Walker->type->Z != 1) // skip hydrogen
2488 #endif
2489 {
2490 AtomMask[Walker->nr] = true; // include all (non-hydrogen) atoms
2491 status = true;
2492 }
2493 }
2494 }
2495 Free((void **)&buffer, "molecule::CheckOrderAtSite: *buffer");
2496 // pick a given number of highest values and set AtomMask
2497 } else { // global increase of Bond Order
2498 while (Walker->next != end) {
2499 Walker = Walker->next;
2500 #ifdef ADDHYDROGEN
2501 if (Walker->type->Z != 1) // skip hydrogen
2502 #endif
2503 {
2504 AtomMask[Walker->nr] = true; // include all (non-hydrogen) atoms
2505 if ((Order != 0) && (Walker->AdaptiveOrder < Order) && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
2506 status = true;
2507 }
2508 }
2509 if ((Order == 0) && (AtomMask[AtomCount] == true)) // single stepping, just check
2510 status = false;
2511
2512 if (!status)
2513 *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
2514 }
2515
2516 // print atom mask for debugging
2517 *out << " ";
2518 for(int i=0;i<AtomCount;i++)
2519 *out << (i % 10);
2520 *out << endl << "Atom mask is: ";
2521 for(int i=0;i<AtomCount;i++)
2522 *out << (AtomMask[i] ? "t" : "f");
2523 *out << endl;
2524
2525 return status;
2526};
2527
2528/** Create a SortIndex to map from BFS labels to the sequence in which the atoms are given in the config file.
2529 * \param *out output stream for debugging
2530 * \param *&SortIndex Mapping array of size molecule::AtomCount
2531 * \return true - success, false - failure of SortIndex alloc
2532 */
2533bool molecule::CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex)
2534{
2535 element *runner = elemente->start;
2536 int AtomNo = 0;
2537 atom *Walker = NULL;
2538
2539 if (SortIndex != NULL) {
2540 *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
2541 return false;
2542 }
2543 SortIndex = (int *) Malloc(sizeof(int)*AtomCount, "molecule::FragmentMolecule: *SortIndex");
2544 for(int i=AtomCount;i--;)
2545 SortIndex[i] = -1;
2546 while (runner->next != elemente->end) { // go through every element
2547 runner = runner->next;
2548 if (ElementsInMolecule[runner->Z]) { // if this element got atoms
2549 Walker = start;
2550 while (Walker->next != end) { // go through every atom of this element
2551 Walker = Walker->next;
2552 if (Walker->type->Z == runner->Z) // if this atom fits to element
2553 SortIndex[Walker->nr] = AtomNo++;
2554 }
2555 }
2556 }
2557 return true;
2558};
2559
2560/** Performs a many-body bond order analysis for a given bond order.
2561 * -# parses adjacency, keysets and orderatsite files
2562 * -# performs DFS to find connected subgraphs (to leave this in was a design decision: might be useful later)
2563 * -# RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energ
2564y contribution", and that's why this consciously not done in the following loop)
2565 * -# in a loop over all subgraphs
2566 * -# calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
2567 * -# creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)
2568 * -# combines the generated molecule lists from all subgraphs
2569 * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
2570 * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
2571 * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
2572 * subgraph in the MoleculeListClass.
2573 * \param *out output stream for debugging
2574 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
2575 * \param *configuration configuration for writing config files for each fragment
2576 */
2577void molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
2578{
2579 MoleculeListClass *BondFragments = NULL;
2580 int *SortIndex = NULL;
2581 int *MinimumRingSize = NULL;
2582 int FragmentCounter;
2583 MoleculeLeafClass *MolecularWalker = NULL;
2584 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
2585 fstream File;
2586 bool FragmentationToDo = true;
2587 Graph **FragmentList = NULL;
2588 Graph *ParsedFragmentList = NULL;
2589 Graph TotalGraph; // graph with all keysets however local numbers
2590 int TotalNumberOfKeySets = 0;
2591 atom **ListOfAtoms = NULL;
2592 atom ***ListOfLocalAtoms = NULL;
2593 bool *AtomMask = NULL;
2594
2595 *out << endl;
2596#ifdef ADDHYDROGEN
2597 *out << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl;
2598#else
2599 *out << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl;
2600#endif
2601
2602 // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
2603
2604 // ===== 1. Check whether bond structure is same as stored in files ====
2605
2606 // fill the adjacency list
2607 CreateListOfBondsPerAtom(out);
2608
2609 // create lookup table for Atom::nr
2610 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);
2611
2612 // === compare it with adjacency file ===
2613 FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
2614 Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
2615
2616 // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
2617 Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
2618 // fill the bond structure of the individually stored subgraphs
2619 Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms
2620
2621 // ===== 3. if structure still valid, parse key set file and others =====
2622 FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
2623
2624 // ===== 4. check globally whether there's something to do actually (first adaptivity check)
2625 FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
2626
2627 // =================================== Begin of FRAGMENTATION ===============================
2628 // ===== 6a. assign each keyset to its respective subgraph =====
2629 Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), false);
2630
2631 KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
2632 AtomMask = new bool[AtomCount+1];
2633 while (CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath)) {
2634 AtomMask[AtomCount] = true; // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
2635 // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
2636 Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
2637
2638 // ===== 7. fill the bond fragment list =====
2639 FragmentCounter = 0;
2640 MolecularWalker = Subgraphs;
2641 while (MolecularWalker->next != NULL) {
2642 MolecularWalker = MolecularWalker->next;
2643 *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
2644 // output ListOfBondsPerAtom for debugging
2645 MolecularWalker->Leaf->OutputListOfBonds(out);
2646 if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
2647
2648 // call BOSSANOVA method
2649 *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
2650 MolecularWalker->Leaf->FragmentBOSSANOVA(out, FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
2651 } else {
2652 cerr << "Subgraph " << MolecularWalker << " has no atoms!" << endl;
2653 }
2654 FragmentCounter++; // next fragment list
2655 }
2656 }
2657 delete[](RootStack);
2658 delete[](AtomMask);
2659 delete(ParsedFragmentList);
2660 delete[](MinimumRingSize);
2661
2662 // free the index lookup list
2663 for (int i=FragmentCounter;i--;)
2664 Free((void **)&ListOfLocalAtoms[i], "molecule::FragmentMolecule - *ListOfLocalAtoms[]");
2665 Free((void **)&ListOfLocalAtoms, "molecule::FragmentMolecule - **ListOfLocalAtoms");
2666
2667 // ==================================== End of FRAGMENTATION ============================================
2668
2669 // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
2670 Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
2671
2672 // free subgraph memory again
2673 FragmentCounter = 0;
2674 if (Subgraphs != NULL) {
2675 while (Subgraphs->next != NULL) {
2676 Subgraphs = Subgraphs->next;
2677 delete(FragmentList[FragmentCounter++]);
2678 delete(Subgraphs->previous);
2679 }
2680 delete(Subgraphs);
2681 }
2682 Free((void **)&FragmentList, "molecule::FragmentMolecule - **FragmentList");
2683
2684 // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
2685 // allocate memory for the pointer array and transmorph graphs into full molecular fragments
2686 BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
2687 int k=0;
2688 for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
2689 KeySet test = (*runner).first;
2690 *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
2691 BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
2692 k++;
2693 }
2694 *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
2695
2696 // ===== 9. Save fragments' configuration and keyset files et al to disk ===
2697 if (BondFragments->NumberOfMolecules != 0) {
2698 // create the SortIndex from BFS labels to order in the config file
2699 CreateMappingLabelsToConfigSequence(out, SortIndex);
2700
2701 *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
2702 if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
2703 *out << Verbose(1) << "All configs written." << endl;
2704 else
2705 *out << Verbose(1) << "Some config writing failed." << endl;
2706
2707 // store force index reference file
2708 BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
2709
2710 // store keysets file
2711 StoreKeySetFile(out, TotalGraph, configuration->configpath);
2712
2713 // store Adjacency file
2714 StoreAdjacencyToFile(out, configuration->configpath);
2715
2716 // store Hydrogen saturation correction file
2717 BondFragments->AddHydrogenCorrection(out, configuration->configpath);
2718
2719 // store adaptive orders into file
2720 StoreOrderAtSiteFile(out, configuration->configpath);
2721
2722 // restore orbital and Stop values
2723 CalculateOrbitals(*configuration);
2724
2725 // free memory for bond part
2726 *out << Verbose(1) << "Freeing bond memory" << endl;
2727 delete(FragmentList); // remove bond molecule from memory
2728 Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
2729 } else
2730 *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
2731
2732 *out << Verbose(0) << "End of bond fragmentation." << endl;
2733};
2734
2735/** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
2736 * Atoms not present in the file get "-1".
2737 * \param *out output stream for debugging
2738 * \param *path path to file ORDERATSITEFILE
2739 * \return true - file writable, false - not writable
2740 */
2741bool molecule::StoreOrderAtSiteFile(ofstream *out, char *path)
2742{
2743 stringstream line;
2744 ofstream file;
2745
2746 line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
2747 file.open(line.str().c_str());
2748 *out << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
2749 if (file != NULL) {
2750 atom *Walker = start;
2751 while (Walker->next != end) {
2752 Walker = Walker->next;
2753 file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << endl;
2754 *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "." << endl;
2755 }
2756 file.close();
2757 *out << Verbose(1) << "done." << endl;
2758 return true;
2759 } else {
2760 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
2761 return false;
2762 }
2763};
2764
2765/** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
2766 * Atoms not present in the file get "0".
2767 * \param *out output stream for debugging
2768 * \param *path path to file ORDERATSITEFILEe
2769 * \return true - file found and scanned, false - file not found
2770 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
2771 */
2772bool molecule::ParseOrderAtSiteFromFile(ofstream *out, char *path)
2773{
2774 unsigned char *OrderArray = (unsigned char *) Malloc(sizeof(unsigned char)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
2775 bool status;
2776 int AtomNr;
2777 stringstream line;
2778 ifstream file;
2779 int Order;
2780
2781 *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
2782 for(int i=AtomCount;i--;)
2783 OrderArray[i] = 0;
2784 line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
2785 file.open(line.str().c_str());
2786 if (file != NULL) {
2787 for (int i=AtomCount;i--;) // initialise with 0
2788 OrderArray[i] = 0;
2789 while (!file.eof()) { // parse from file
2790 file >> AtomNr;
2791 file >> Order;
2792 OrderArray[AtomNr] = (unsigned char) Order;
2793 //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << "." << endl;
2794 }
2795 atom *Walker = start;
2796 while (Walker->next != end) { // fill into atom classes
2797 Walker = Walker->next;
2798 Walker->AdaptiveOrder = OrderArray[Walker->nr];
2799 *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << "." << endl;
2800 }
2801 file.close();
2802 *out << Verbose(1) << "done." << endl;
2803 status = true;
2804 } else {
2805 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
2806 status = false;
2807 }
2808 Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
2809
2810 *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
2811 return status;
2812};
2813
2814/** Creates an 2d array of pointer with an entry for each atom and each bond it has.
2815 * Updates molecule::ListOfBondsPerAtom, molecule::NumberOfBondsPerAtom by parsing through
2816 * bond chain list, using molecule::AtomCount and molecule::BondCount.
2817 * Allocates memory, fills the array and exits
2818 * \param *out output stream for debugging
2819 */
2820void molecule::CreateListOfBondsPerAtom(ofstream *out)
2821{
2822 bond *Binder = NULL;
2823 atom *Walker = NULL;
2824 int TotalDegree;
2825 *out << Verbose(1) << "Begin of Creating ListOfBondsPerAtom: AtomCount = " << AtomCount << "\tBondCount = " << BondCount << "\tNoNonBonds = " << NoNonBonds << "." << endl;
2826
2827 // re-allocate memory
2828 *out << Verbose(2) << "(Re-)Allocating memory." << endl;
2829 if (ListOfBondsPerAtom != NULL) {
2830 for(int i=AtomCount;i--;)
2831 Free((void **)&ListOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: ListOfBondsPerAtom[i]");
2832 Free((void **)&ListOfBondsPerAtom, "molecule::CreateListOfBondsPerAtom: ListOfBondsPerAtom");
2833 }
2834 if (NumberOfBondsPerAtom != NULL)
2835 Free((void **)&NumberOfBondsPerAtom, "molecule::CreateListOfBondsPerAtom: NumberOfBondsPerAtom");
2836 ListOfBondsPerAtom = (bond ***) Malloc(sizeof(bond **)*AtomCount, "molecule::CreateListOfBondsPerAtom: ***ListOfBondsPerAtom");
2837 NumberOfBondsPerAtom = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfBondsPerAtom: *NumberOfBondsPerAtom");
2838
2839 // reset bond counts per atom
2840 for(int i=AtomCount;i--;)
2841 NumberOfBondsPerAtom[i] = 0;
2842 // count bonds per atom
2843 Binder = first;
2844 while (Binder->next != last) {
2845 Binder = Binder->next;
2846 NumberOfBondsPerAtom[Binder->leftatom->nr]++;
2847 NumberOfBondsPerAtom[Binder->rightatom->nr]++;
2848 }
2849 for(int i=AtomCount;i--;) {
2850 // allocate list of bonds per atom
2851 ListOfBondsPerAtom[i] = (bond **) Malloc(sizeof(bond *)*NumberOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: **ListOfBondsPerAtom[]");
2852 // clear the list again, now each NumberOfBondsPerAtom marks current free field
2853 NumberOfBondsPerAtom[i] = 0;
2854 }
2855 // fill the list
2856 Binder = first;
2857 while (Binder->next != last) {
2858 Binder = Binder->next;
2859 ListOfBondsPerAtom[Binder->leftatom->nr][NumberOfBondsPerAtom[Binder->leftatom->nr]++] = Binder;
2860 ListOfBondsPerAtom[Binder->rightatom->nr][NumberOfBondsPerAtom[Binder->rightatom->nr]++] = Binder;
2861 }
2862
2863 // output list for debugging
2864 *out << Verbose(3) << "ListOfBondsPerAtom for each atom:" << endl;
2865 Walker = start;
2866 while (Walker->next != end) {
2867 Walker = Walker->next;
2868 *out << Verbose(4) << "Atom " << Walker->Name << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
2869 TotalDegree = 0;
2870 for (int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {
2871 *out << *ListOfBondsPerAtom[Walker->nr][j] << "\t";
2872 TotalDegree += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
2873 }
2874 *out << " -- TotalDegree: " << TotalDegree << endl;
2875 }
2876 *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
2877};
2878
2879/** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
2880 * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
2881 * white and putting into queue.
2882 * \param *out output stream for debugging
2883 * \param *Mol Molecule class to add atoms to
2884 * \param **AddedAtomList list with added atom pointers, index is atom father's number
2885 * \param **AddedBondList list with added bond pointers, index is bond father's number
2886 * \param *Root root vertex for BFS
2887 * \param *Bond bond not to look beyond
2888 * \param BondOrder maximum distance for vertices to add
2889 * \param IsAngstroem lengths are in angstroem or bohrradii
2890 */
2891void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
2892{
2893 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::BreadthFirstSearchAdd: **PredecessorList");
2894 int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::BreadthFirstSearchAdd: *ShortestPathList");
2895 enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::BreadthFirstSearchAdd: *ColorList");
2896 class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
2897 atom *Walker = NULL, *OtherAtom = NULL;
2898 bond *Binder = NULL;
2899
2900 // add Root if not done yet
2901 AtomStack->ClearStack();
2902 if (AddedAtomList[Root->nr] == NULL) // add Root if not yet present
2903 AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
2904 AtomStack->Push(Root);
2905
2906 // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
2907 for (int i=AtomCount;i--;) {
2908 PredecessorList[i] = NULL;
2909 ShortestPathList[i] = -1;
2910 if (AddedAtomList[i] != NULL) // mark already present atoms (i.e. Root and maybe others) as visited
2911 ColorList[i] = lightgray;
2912 else
2913 ColorList[i] = white;
2914 }
2915 ShortestPathList[Root->nr] = 0;
2916
2917 // and go on ... Queue always contains all lightgray vertices
2918 while (!AtomStack->IsEmpty()) {
2919 // we have to pop the oldest atom from stack. This keeps the atoms on the stack always of the same ShortestPath distance.
2920 // 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
2921 // append length of 3 (their neighbours). Thus on stack we have always atoms of a certain length n at bottom of stack and
2922 // followed by n+1 till top of stack.
2923 Walker = AtomStack->PopFirst(); // pop oldest added
2924 *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
2925 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
2926 Binder = ListOfBondsPerAtom[Walker->nr][i];
2927 if (Binder != NULL) { // don't look at bond equal NULL
2928 OtherAtom = Binder->GetOtherAtom(Walker);
2929 *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
2930 if (ColorList[OtherAtom->nr] == white) {
2931 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)
2932 ColorList[OtherAtom->nr] = lightgray;
2933 PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
2934 ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
2935 *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;
2936 if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
2937 *out << Verbose(3);
2938 if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
2939 AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
2940 *out << "Added OtherAtom " << OtherAtom->Name;
2941 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2942 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2943 AddedBondList[Binder->nr]->Type = Binder->Type;
2944 *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
2945 } 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)
2946 *out << "Not adding OtherAtom " << OtherAtom->Name;
2947 if (AddedBondList[Binder->nr] == NULL) {
2948 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2949 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2950 AddedBondList[Binder->nr]->Type = Binder->Type;
2951 *out << ", added Bond " << *(AddedBondList[Binder->nr]);
2952 } else
2953 *out << ", not added Bond ";
2954 }
2955 *out << ", putting OtherAtom into queue." << endl;
2956 AtomStack->Push(OtherAtom);
2957 } else { // out of bond order, then replace
2958 if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
2959 ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
2960 if (Binder == Bond)
2961 *out << Verbose(3) << "Not Queueing, is the Root bond";
2962 else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
2963 *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
2964 if (!Binder->Cyclic)
2965 *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
2966 if (AddedBondList[Binder->nr] == NULL) {
2967 if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
2968 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2969 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2970 AddedBondList[Binder->nr]->Type = Binder->Type;
2971 } else {
2972#ifdef ADDHYDROGEN
2973 Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem);
2974#endif
2975 }
2976 }
2977 }
2978 } else {
2979 *out << Verbose(3) << "Not Adding, has already been visited." << endl;
2980 // This has to be a cyclic bond, check whether it's present ...
2981 if (AddedBondList[Binder->nr] == NULL) {
2982 if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
2983 AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
2984 AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
2985 AddedBondList[Binder->nr]->Type = Binder->Type;
2986 } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
2987#ifdef ADDHYDROGEN
2988 Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem);
2989#endif
2990 }
2991 }
2992 }
2993 }
2994 }
2995 ColorList[Walker->nr] = black;
2996 *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
2997 }
2998 Free((void **)&PredecessorList, "molecule::BreadthFirstSearchAdd: **PredecessorList");
2999 Free((void **)&ShortestPathList, "molecule::BreadthFirstSearchAdd: **ShortestPathList");
3000 Free((void **)&ColorList, "molecule::BreadthFirstSearchAdd: **ColorList");
3001 delete(AtomStack);
3002};
3003
3004/** Adds bond structure to this molecule from \a Father molecule.
3005 * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
3006 * with end points present in this molecule, bond is created in this molecule.
3007 * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
3008 * \param *out output stream for debugging
3009 * \param *Father father molecule
3010 * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
3011 * \todo not checked, not fully working probably
3012 */
3013bool molecule::BuildInducedSubgraph(ofstream *out, const molecule *Father)
3014{
3015 atom *Walker = NULL, *OtherAtom = NULL;
3016 bool status = true;
3017 atom **ParentList = (atom **) Malloc(sizeof(atom *)*Father->AtomCount, "molecule::BuildInducedSubgraph: **ParentList");
3018
3019 *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
3020
3021 // reset parent list
3022 *out << Verbose(3) << "Resetting ParentList." << endl;
3023 for (int i=Father->AtomCount;i--;)
3024 ParentList[i] = NULL;
3025
3026 // fill parent list with sons
3027 *out << Verbose(3) << "Filling Parent List." << endl;
3028 Walker = start;
3029 while (Walker->next != end) {
3030 Walker = Walker->next;
3031 ParentList[Walker->father->nr] = Walker;
3032 // Outputting List for debugging
3033 *out << Verbose(4) << "Son["<< Walker->father->nr <<"] of " << Walker->father << " is " << ParentList[Walker->father->nr] << "." << endl;
3034 }
3035
3036 // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
3037 *out << Verbose(3) << "Creating bonds." << endl;
3038 Walker = Father->start;
3039 while (Walker->next != Father->end) {
3040 Walker = Walker->next;
3041 if (ParentList[Walker->nr] != NULL) {
3042 if (ParentList[Walker->nr]->father != Walker) {
3043 status = false;
3044 } else {
3045 for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
3046 OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
3047 if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
3048 *out << Verbose(4) << "Endpoints of Bond " << Father->ListOfBondsPerAtom[Walker->nr][i] << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
3049 AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], Father->ListOfBondsPerAtom[Walker->nr][i]->BondDegree);
3050 }
3051 }
3052 }
3053 }
3054 }
3055
3056 Free((void **)&ParentList, "molecule::BuildInducedSubgraph: **ParentList");
3057 *out << Verbose(2) << "End of BuildInducedSubgraph." << endl;
3058 return status;
3059};
3060
3061
3062/** Looks through a StackClass<atom *> and returns the likeliest removal candiate.
3063 * \param *out output stream for debugging messages
3064 * \param *&Leaf KeySet to look through
3065 * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
3066 * \param index of the atom suggested for removal
3067 */
3068int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
3069{
3070 atom *Runner = NULL;
3071 int SP, Removal;
3072
3073 *out << Verbose(2) << "Looking for removal candidate." << endl;
3074 SP = -1; //0; // not -1, so that Root is never removed
3075 Removal = -1;
3076 for (KeySet::iterator runner = Leaf->begin(); runner != Leaf->end(); runner++) {
3077 Runner = FindAtom((*runner));
3078 if (Runner->type->Z != 1) { // skip all those added hydrogens when re-filling snake stack
3079 if (ShortestPathList[(*runner)] > SP) { // remove the oldest one with longest shortest path
3080 SP = ShortestPathList[(*runner)];
3081 Removal = (*runner);
3082 }
3083 }
3084 }
3085 return Removal;
3086};
3087
3088/** Stores a fragment from \a KeySet into \a molecule.
3089 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
3090 * molecule and adds missing hydrogen where bonds were cut.
3091 * \param *out output stream for debugging messages
3092 * \param &Leaflet pointer to KeySet structure
3093 * \param IsAngstroem whether we have Ansgtroem or bohrradius
3094 * \return pointer to constructed molecule
3095 */
3096molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
3097{
3098 atom *Runner = NULL, *FatherOfRunner = NULL, *OtherFather = NULL;
3099 atom **SonList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::StoreFragmentFromStack: **SonList");
3100 molecule *Leaf = new molecule(elemente);
3101
3102// *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
3103
3104 Leaf->BondDistance = BondDistance;
3105 for(int i=NDIM*2;i--;)
3106 Leaf->cell_size[i] = cell_size[i];
3107
3108 // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
3109 for(int i=AtomCount;i--;)
3110 SonList[i] = NULL;
3111
3112 // first create the minimal set of atoms from the KeySet
3113 for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
3114 FatherOfRunner = FindAtom((*runner)); // find the id
3115 SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
3116 }
3117
3118 // create the bonds between all: Make it an induced subgraph and add hydrogen
3119// *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
3120 Runner = Leaf->start;
3121 while (Runner->next != Leaf->end) {
3122 Runner = Runner->next;
3123 FatherOfRunner = Runner->father;
3124 if (SonList[FatherOfRunner->nr] != NULL) { // check if this, our father, is present in list
3125 // create all bonds
3126 for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
3127 OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
3128// *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
3129 if (SonList[OtherFather->nr] != NULL) {
3130// *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
3131 if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
3132// *out << Verbose(3) << "Adding Bond: ";
3133// *out <<
3134 Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
3135// *out << "." << endl;
3136 //NumBonds[Runner->nr]++;
3137 } else {
3138// *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
3139 }
3140 } else {
3141// *out << ", who has no son in this fragment molecule." << endl;
3142#ifdef ADDHYDROGEN
3143 //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
3144 Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem);
3145#endif
3146 //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
3147 }
3148 }
3149 } else {
3150 *out << Verbose(0) << "ERROR: Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
3151 }
3152#ifdef ADDHYDROGEN
3153 while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
3154 Runner = Runner->next;
3155#endif
3156 }
3157 Leaf->CreateListOfBondsPerAtom(out);
3158 //Leaflet->Leaf->ScanForPeriodicCorrection(out);
3159 Free((void **)&SonList, "molecule::StoreFragmentFromStack: **SonList");
3160// *out << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
3161 return Leaf;
3162};
3163
3164/** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
3165 * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
3166 * computer game, that winds through the connected graph representing the molecule. Color (white,
3167 * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
3168 * creating only unique fragments and not additional ones with vertices simply in different sequence.
3169 * The Predecessor is always the one that came before in discovering, needed on backstepping. And
3170 * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
3171 * stepping.
3172 * \param *out output stream for debugging
3173 * \param Order number of atoms in each fragment
3174 * \param *configuration configuration for writing config files for each fragment
3175 * \return List of all unique fragments with \a Order atoms
3176 */
3177/*
3178MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
3179{
3180 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
3181 int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
3182 int *Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
3183 enum Shading *ColorVertexList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
3184 enum Shading *ColorEdgeList = (enum Shading *) Malloc(sizeof(enum Shading)*BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
3185 StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
3186 StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
3187 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!
3188 MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
3189 MoleculeListClass *FragmentList = NULL;
3190 atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
3191 bond *Binder = NULL;
3192 int RunningIndex = 0, FragmentCounter = 0;
3193
3194 *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
3195
3196 // reset parent list
3197 *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
3198 for (int i=0;i<AtomCount;i++) { // reset all atom labels
3199 // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
3200 Labels[i] = -1;
3201 SonList[i] = NULL;
3202 PredecessorList[i] = NULL;
3203 ColorVertexList[i] = white;
3204 ShortestPathList[i] = -1;
3205 }
3206 for (int i=0;i<BondCount;i++)
3207 ColorEdgeList[i] = white;
3208 RootStack->ClearStack(); // clearstack and push first atom if exists
3209 TouchedStack->ClearStack();
3210 Walker = start->next;
3211 while ((Walker != end)
3212#ifdef ADDHYDROGEN
3213 && (Walker->type->Z == 1)
3214#endif
3215 ) { // search for first non-hydrogen atom
3216 *out << Verbose(4) << "Current Root candidate is " << Walker->Name << "." << endl;
3217 Walker = Walker->next;
3218 }
3219 if (Walker != end)
3220 RootStack->Push(Walker);
3221 else
3222 *out << Verbose(0) << "ERROR: Could not find an appropriate Root atom!" << endl;
3223 *out << Verbose(3) << "Root " << Walker->Name << " is on AtomStack, beginning loop through all vertices ..." << endl;
3224
3225 ///// OUTER LOOP ////////////
3226 while (!RootStack->IsEmpty()) {
3227 // get new root vertex from atom stack
3228 Root = RootStack->PopFirst();
3229 ShortestPathList[Root->nr] = 0;
3230 if (Labels[Root->nr] == -1)
3231 Labels[Root->nr] = RunningIndex++; // prevent it from getting again on AtomStack
3232 PredecessorList[Root->nr] = Root;
3233 TouchedStack->Push(Root);
3234 *out << Verbose(0) << "Root for this loop is: " << Root->Name << ".\n";
3235
3236 // clear snake stack
3237 SnakeStack->ClearStack();
3238 //SnakeStack->TestImplementation(out, start->next);
3239
3240 ///// INNER LOOP ////////////
3241 // Problems:
3242 // - what about cyclic bonds?
3243 Walker = Root;
3244 do {
3245 *out << Verbose(1) << "Current Walker is: " << Walker->Name;
3246 // initial setting of the new Walker: label, color, shortest path and put on stacks
3247 if (Labels[Walker->nr] == -1) { // give atom a unique, monotonely increasing number
3248 Labels[Walker->nr] = RunningIndex++;
3249 RootStack->Push(Walker);
3250 }
3251 *out << ", has label " << Labels[Walker->nr];
3252 if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) { // color it if newly discovered and push on stacks (and if within reach!)
3253 if ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white)) {
3254 // Binder ought to be set still from last neighbour search
3255 *out << ", coloring bond " << *Binder << " black";
3256 ColorEdgeList[Binder->nr] = black; // mark this bond as used
3257 }
3258 if (ShortestPathList[Walker->nr] == -1) {
3259 ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1;
3260 TouchedStack->Push(Walker); // mark every atom for lists cleanup later, whose shortest path has been changed
3261 }
3262 if ((ShortestPathList[Walker->nr] < Order) && (ColorVertexList[Walker->nr] != darkgray)) { // if not already on snake stack
3263 SnakeStack->Push(Walker);
3264 ColorVertexList[Walker->nr] = darkgray; // mark as dark gray of on snake stack
3265 }
3266 }
3267 *out << ", SP of " << ShortestPathList[Walker->nr] << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
3268
3269 // then check the stack for a newly stumbled upon fragment
3270 if (SnakeStack->ItemCount() == Order) { // is stack full?
3271 // store the fragment if it is one and get a removal candidate
3272 Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
3273 // remove the candidate if one was found
3274 if (Removal != NULL) {
3275 *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
3276 SnakeStack->RemoveItem(Removal);
3277 ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
3278 if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
3279 Walker = PredecessorList[Removal->nr];
3280 *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
3281 }
3282 }
3283 } else
3284 Removal = NULL;
3285
3286 // finally, look for a white neighbour as the next Walker
3287 Binder = NULL;
3288 if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) { // don't look, if a new walker has been set above
3289 *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
3290 OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
3291 if (ShortestPathList[Walker->nr] < Order) {
3292 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
3293 Binder = ListOfBondsPerAtom[Walker->nr][i];
3294 *out << Verbose(2) << "Current bond is " << *Binder << ": ";
3295 OtherAtom = Binder->GetOtherAtom(Walker);
3296 if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
3297 *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
3298 //ColorVertexList[OtherAtom->nr] = lightgray; // mark as explored
3299 } else { // otherwise check its colour and element
3300 if (
3301#ifdef ADDHYDROGEN
3302 (OtherAtom->type->Z != 1) &&
3303#endif
3304 (ColorEdgeList[Binder->nr] == white)) { // skip hydrogen, look for unexplored vertices
3305 *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
3306 // i find it currently rather sensible to always set the predecessor in order to find one's way back
3307 //if (PredecessorList[OtherAtom->nr] == NULL) {
3308 PredecessorList[OtherAtom->nr] = Walker;
3309 *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
3310 //} else {
3311 // *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
3312 //}
3313 Walker = OtherAtom;
3314 break;
3315 } else {
3316 if (OtherAtom->type->Z == 1)
3317 *out << "Links to a hydrogen atom." << endl;
3318 else
3319 *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
3320 }
3321 }
3322 }
3323 } else { // means we have stepped beyond the horizon: Return!
3324 Walker = PredecessorList[Walker->nr];
3325 OtherAtom = Walker;
3326 *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
3327 }
3328 if (Walker != OtherAtom) { // if no white neighbours anymore, color it black
3329 *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
3330 ColorVertexList[Walker->nr] = black;
3331 Walker = PredecessorList[Walker->nr];
3332 }
3333 }
3334 } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
3335 *out << Verbose(2) << "Inner Looping is finished." << endl;
3336
3337 // if we reset all AtomCount atoms, we have again technically O(N^2) ...
3338 *out << Verbose(2) << "Resetting lists." << endl;
3339 Walker = NULL;
3340 Binder = NULL;
3341 while (!TouchedStack->IsEmpty()) {
3342 Walker = TouchedStack->PopLast();
3343 *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
3344 for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
3345 ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
3346 PredecessorList[Walker->nr] = NULL;
3347 ColorVertexList[Walker->nr] = white;
3348 ShortestPathList[Walker->nr] = -1;
3349 }
3350 }
3351 *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
3352
3353 // copy together
3354 *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
3355 FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
3356 RunningIndex = 0;
3357 while ((Leaflet != NULL) && (RunningIndex < FragmentCounter)) {
3358 FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
3359 Leaflet->Leaf = NULL; // prevent molecule from being removed
3360 TempLeaf = Leaflet;
3361 Leaflet = Leaflet->previous;
3362 delete(TempLeaf);
3363 };
3364
3365 // free memory and exit
3366 Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
3367 Free((void **)&ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
3368 Free((void **)&Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
3369 Free((void **)&ColorVertexList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
3370 delete(RootStack);
3371 delete(TouchedStack);
3372 delete(SnakeStack);
3373
3374 *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
3375 return FragmentList;
3376};
3377*/
3378
3379/** Structure containing all values in power set combination generation.
3380 */
3381struct UniqueFragments {
3382 config *configuration;
3383 atom *Root;
3384 Graph *Leaflet;
3385 KeySet *FragmentSet;
3386 int ANOVAOrder;
3387 int FragmentCounter;
3388 int CurrentIndex;
3389 int *Labels;
3390 double TEFactor;
3391 int *ShortestPathList;
3392 bool **UsedList;
3393 bond **BondsPerSPList;
3394 int *BondsPerSPCount;
3395};
3396
3397/** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.
3398 * -# loops over every possible combination (2^dimension of edge set)
3399 * -# inserts current set, if there's still space left
3400 * -# yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root dist
3401ance+1
3402 * -# no: stores fragment into keyset list by calling InsertFragmentIntoGraph
3403 * -# removes all items added into the snake stack (in UniqueFragments structure) added during level (root
3404distance) and current set
3405 * \param *out output stream for debugging
3406 * \param FragmentSearch UniqueFragments structure with all values needed
3407 * \param RootDistance current shortest path level, whose set of edges is represented by **BondsSet
3408 * \param SetDimension Number of possible bonds on this level (i.e. size of the array BondsSet[])
3409 * \param SubOrder remaining number of allowed vertices to add
3410 */
3411void molecule::SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder)
3412{
3413 atom *OtherWalker = NULL;
3414 int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
3415 int NumCombinations;
3416 bool bit;
3417 int bits, TouchedIndex, SubSetDimension, SP;
3418 int Removal;
3419 int *TouchedList = (int *) Malloc(sizeof(int)*(SubOrder+1), "molecule::SPFragmentGenerator: *TouchedList");
3420 bond *Binder = NULL;
3421 bond **BondsList = NULL;
3422
3423 NumCombinations = 1 << SetDimension;
3424
3425 // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
3426 // von Endstuecken (aus den Bonds) hinzugefᅵᅵgt werden und fᅵᅵr verbleibende ANOVAOrder
3427 // rekursiv GraphCrawler in der nᅵᅵchsten Ebene aufgerufen werden
3428
3429 *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
3430 *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;
3431
3432 // initialised touched list (stores added atoms on this level)
3433 *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
3434 for (TouchedIndex=SubOrder+1;TouchedIndex--;) // empty touched list
3435 TouchedList[TouchedIndex] = -1;
3436 TouchedIndex = 0;
3437
3438 // create every possible combination of the endpieces
3439 *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
3440 for (int i=1;i<NumCombinations;i++) { // sweep through all power set combinations (skip empty set!)
3441 // count the set bit of i
3442 bits = 0;
3443 for (int j=SetDimension;j--;)
3444 bits += (i & (1 << j)) >> j;
3445
3446 *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
3447 if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
3448 // --1-- add this set of the power set of bond partners to the snake stack
3449 for (int j=0;j<SetDimension;j++) { // pull out every bit by shifting
3450 bit = ((i & (1 << j)) != 0); // mask the bit for the j-th bond
3451 if (bit) { // if bit is set, we add this bond partner
3452 OtherWalker = BondsSet[j]->rightatom; // rightatom is always the one more distant, i.e. the one to add
3453 //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
3454 //if ((!FragmentSearch->UsedList[OtherWalker->nr][i]) && (FragmentSearch->Labels[OtherWalker->nr] > FragmentSearch->Labels[FragmentSearch->Root->nr])) {
3455 //*out << Verbose(2+verbosity) << "Not used so far, label " << FragmentSearch->Labels[OtherWalker->nr] << " is bigger than Root's " << FragmentSearch->Labels[FragmentSearch->Root->nr] << "." << endl;
3456 *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
3457 TouchedList[TouchedIndex++] = OtherWalker->nr; // note as added
3458 FragmentSearch->FragmentSet->insert(OtherWalker->nr);
3459 //FragmentSearch->UsedList[OtherWalker->nr][i] = true;
3460 //}
3461 } else {
3462 *out << Verbose(2+verbosity) << "Not adding." << endl;
3463 }
3464 }
3465
3466 if (bits < SubOrder) {
3467 *out << Verbose(1+verbosity) << "There's still some space left on stack: " << (SubOrder - bits) << "." << endl;
3468 // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
3469 SP = RootDistance+1; // this is the next level
3470 // first count the members in the subset
3471 SubSetDimension = 0;
3472 Binder = FragmentSearch->BondsPerSPList[2*SP]; // start node for this level
3473 while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) { // compare to end node of this level
3474 Binder = Binder->next;
3475 for (int k=TouchedIndex;k--;) {
3476 if (Binder->Contains(TouchedList[k])) // if we added this very endpiece
3477 SubSetDimension++;
3478 }
3479 }
3480 // then allocate and fill the list
3481 BondsList = (bond **) Malloc(sizeof(bond *)*SubSetDimension, "molecule::SPFragmentGenerator: **BondsList");
3482 SubSetDimension = 0;
3483 Binder = FragmentSearch->BondsPerSPList[2*SP];
3484 while (Binder->next != FragmentSearch->BondsPerSPList[2*SP+1]) {
3485 Binder = Binder->next;
3486 for (int k=0;k<TouchedIndex;k++) {
3487 if (Binder->leftatom->nr == TouchedList[k]) // leftatom is always the close one
3488 BondsList[SubSetDimension++] = Binder;
3489 }
3490 }
3491 *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
3492 SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
3493 Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
3494 } else {
3495 // --2-- otherwise store the complete fragment
3496 *out << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
3497 // store fragment as a KeySet
3498 *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
3499 for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
3500 *out << (*runner) << " ";
3501 InsertFragmentIntoGraph(out, FragmentSearch);
3502 //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
3503 //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList,FragmentSearch->Labels, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
3504 }
3505
3506 // --3-- remove all added items in this level from snake stack
3507 *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
3508 for(int j=0;j<TouchedIndex;j++) {
3509 Removal = TouchedList[j];
3510 *out << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
3511 FragmentSearch->FragmentSet->erase(Removal);
3512 TouchedList[j] = -1;
3513 }
3514 *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
3515 for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
3516 *out << (*runner) << " ";
3517 *out << endl;
3518 TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
3519 } else {
3520 *out << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
3521 }
3522 }
3523 Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");
3524 *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
3525};
3526
3527/** 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.
3528 * -# initialises UniqueFragments structure
3529 * -# fills edge list via BFS
3530 * -# creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as
3531 root distance, the edge set, its dimension and the current suborder
3532 * -# Free'ing structure
3533 * 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
3534 * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
3535 * \param *out output stream for debugging
3536 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
3537 * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
3538 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
3539 * \return number of inserted fragments
3540 * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
3541 */
3542int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
3543{
3544 int SP, UniqueIndex, AtomKeyNr;
3545 int *NumberOfAtomsSPLevel = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *SPLevelCount");
3546 atom *Walker = NULL, *OtherWalker = NULL;
3547 bond *Binder = NULL;
3548 bond **BondsList = NULL;
3549 KeyStack AtomStack;
3550 atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::PowerSetGenerator: **PredecessorList");
3551 int RootKeyNr = FragmentSearch.Root->nr;
3552 int Counter = FragmentSearch.FragmentCounter;
3553
3554 for (int i=AtomCount;i--;) {
3555 PredecessorList[i] = NULL;
3556 }
3557
3558 *out << endl;
3559 *out << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl;
3560
3561 UniqueIndex = 0;
3562 if (FragmentSearch.Labels[RootKeyNr] == -1)
3563 FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
3564 FragmentSearch.ShortestPathList[RootKeyNr] = 0;
3565 // prepare the atom stack counters (number of atoms with certain SP on stack)
3566 for (int i=Order;i--;)
3567 NumberOfAtomsSPLevel[i] = 0;
3568 NumberOfAtomsSPLevel[0] = 1; // for root
3569 SP = -1;
3570 *out << endl;
3571 *out << Verbose(0) << "Starting BFS analysis ..." << endl;
3572 // push as first on atom stack and goooo ...
3573 AtomStack.clear();
3574 AtomStack.push_back(RootKeyNr);
3575 *out << Verbose(0) << "Cleared AtomStack and pushed root as first item onto it." << endl;
3576 // do a BFS search to fill the SP lists and label the found vertices
3577 while (!AtomStack.empty()) {
3578 // pop next atom
3579 AtomKeyNr = AtomStack.front();
3580 AtomStack.pop_front();
3581 if (SP != -1)
3582 NumberOfAtomsSPLevel[SP]--;
3583 if ((SP == -1) || (NumberOfAtomsSPLevel[SP] == -1)) {
3584 SP++;
3585 NumberOfAtomsSPLevel[SP]--; // carry over "-1" to next level
3586 *out << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with 0 item(s)";
3587 if (SP > 0)
3588 *out << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl;
3589 else
3590 *out << "." << endl;
3591 FragmentSearch.BondsPerSPCount[SP] = 0;
3592 } else {
3593 *out << Verbose(1) << "Still " << NumberOfAtomsSPLevel[SP]+1 << " on this SP (" << SP << ") level, currently having " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
3594 }
3595 Walker = FindAtom(AtomKeyNr);
3596 *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << "." << endl;
3597 // check for new sp level
3598 // go through all its bonds
3599 *out << Verbose(1) << "Going through all bonds of Walker." << endl;
3600 for (int i=0;i<NumberOfBondsPerAtom[AtomKeyNr];i++) {
3601 Binder = ListOfBondsPerAtom[AtomKeyNr][i];
3602 OtherWalker = Binder->GetOtherAtom(Walker);
3603 if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
3604#ifdef ADDHYDROGEN
3605 && (OtherWalker->type->Z != 1)
3606#endif
3607 ) { // skip hydrogens and restrict to fragment
3608 *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
3609 // set the label if not set (and push on root stack as well)
3610 if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
3611 FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
3612 *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
3613 } else {
3614 *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
3615 }
3616 if ((OtherWalker != PredecessorList[AtomKeyNr]) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
3617 FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
3618 *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
3619 if (FragmentSearch.ShortestPathList[OtherWalker->nr] < Order) { // only pass through those within reach of Order
3620 // push for exploration on stack (only if the SP of OtherWalker is longer than of Walker! (otherwise it has been added already!)
3621 if (FragmentSearch.ShortestPathList[OtherWalker->nr] > SP) {
3622 if (FragmentSearch.ShortestPathList[OtherWalker->nr] < (Order-1)) {
3623 *out << Verbose(3) << "Putting on atom stack for further exploration." << endl;
3624 PredecessorList[OtherWalker->nr] = Walker; // note down the one further up the exploration tree
3625 AtomStack.push_back(OtherWalker->nr);
3626 NumberOfAtomsSPLevel[FragmentSearch.ShortestPathList[OtherWalker->nr]]++;
3627 } else {
3628 *out << Verbose(3) << "Not putting on atom stack, is at end of reach." << endl;
3629 }
3630 // add the bond in between to the SP list
3631 Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
3632 add(Binder, FragmentSearch.BondsPerSPList[2*SP+1]);
3633 FragmentSearch.BondsPerSPCount[SP]++;
3634 *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP] << " item(s)." << endl;
3635 } else *out << Verbose(3) << "Not putting on atom stack, nor adding bond, as " << *OtherWalker << "'s SP is shorter than Walker's." << endl;
3636 } else *out << Verbose(3) << "Not continuing, as " << *OtherWalker << " is out of reach of order " << Order << "." << endl;
3637 } 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;
3638 } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
3639 }
3640 }
3641 // reset predecessor list
3642 for(int i=0;i<Order;i++) {
3643 Binder = FragmentSearch.BondsPerSPList[2*i];
3644 *out << Verbose(1) << "Current SP level is " << i << "." << endl;
3645 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3646 Binder = Binder->next;
3647 PredecessorList[Binder->rightatom->nr] = NULL; // By construction "OtherAtom" is always Bond::rightatom!
3648 }
3649 }
3650 *out << endl;
3651
3652 // outputting all list for debugging
3653 *out << Verbose(0) << "Printing all found lists." << endl;
3654 for(int i=0;i<Order;i++) {
3655 Binder = FragmentSearch.BondsPerSPList[2*i];
3656 *out << Verbose(1) << "Current SP level is " << i << "." << endl;
3657 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3658 Binder = Binder->next;
3659 *out << Verbose(2) << *Binder << endl;
3660 }
3661 }
3662
3663 // creating fragments with the found edge sets (may be done in reverse order, faster)
3664 SP = 0;
3665 for(int i=Order;i--;) { // sum up all found edges
3666 Binder = FragmentSearch.BondsPerSPList[2*i];
3667 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3668 Binder = Binder->next;
3669 SP ++;
3670 }
3671 }
3672 *out << Verbose(0) << "Total number of edges is " << SP << "." << endl;
3673 if (SP >= (Order-1)) {
3674 // start with root (push on fragment stack)
3675 *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
3676 FragmentSearch.FragmentSet->clear();
3677 FragmentSearch.FragmentSet->insert(FragmentSearch.Root->nr);
3678
3679 if (FragmentSearch.FragmentSet->size() == (unsigned int) Order) {
3680 *out << Verbose(0) << "Enough items on stack already for a fragment!" << endl;
3681 // store fragment as a KeySet
3682 *out << Verbose(2) << "Found a new fragment[" << FragmentSearch.FragmentCounter << "], local nr.s are: ";
3683 for(KeySet::iterator runner = FragmentSearch.FragmentSet->begin(); runner != FragmentSearch.FragmentSet->end(); runner++) {
3684 *out << (*runner) << " ";
3685 }
3686 *out << endl;
3687 InsertFragmentIntoGraph(out, &FragmentSearch);
3688 } else {
3689 *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
3690 // prepare the subset and call the generator
3691 BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
3692 Binder = FragmentSearch.BondsPerSPList[0];
3693 for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++) {
3694 Binder = Binder->next;
3695 BondsList[i] = Binder;
3696 }
3697 SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order-1);
3698 Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
3699 }
3700 } else {
3701 *out << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl;
3702 }
3703
3704 // as FragmentSearch structure is used only once, we don't have to clean it anymore
3705 // remove root from stack
3706 *out << Verbose(0) << "Removing root again from stack." << endl;
3707 FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
3708
3709 // free'ing the bonds lists
3710 *out << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl;
3711 for(int i=Order;i--;) {
3712 *out << Verbose(1) << "Current SP level is " << i << ": ";
3713 Binder = FragmentSearch.BondsPerSPList[2*i];
3714 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
3715 Binder = Binder->next;
3716 // *out << "Removing atom " << Binder->leftatom->nr << " and " << Binder->rightatom->nr << "." << endl; // make sure numbers are local
3717 FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
3718 FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
3719 }
3720 // delete added bonds
3721 cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
3722 // also start and end node
3723 *out << "cleaned." << endl;
3724 }
3725
3726 // free allocated memory
3727 Free((void **)&NumberOfAtomsSPLevel, "molecule::PowerSetGenerator: *SPLevelCount");
3728 Free((void **)&PredecessorList, "molecule::PowerSetGenerator: **PredecessorList");
3729
3730 // return list
3731 *out << Verbose(0) << "End of PowerSetGenerator." << endl;
3732 return (FragmentSearch.FragmentCounter - Counter);
3733};
3734
3735/** Corrects the nuclei position if the fragment was created over the cell borders.
3736 * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
3737 * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
3738 * and re-add the bond. Looping on the distance check.
3739 * \param *out ofstream for debugging messages
3740 */
3741void molecule::ScanForPeriodicCorrection(ofstream *out)
3742{
3743 bond *Binder = NULL;
3744 bond *OtherBinder = NULL;
3745 atom *Walker = NULL;
3746 atom *OtherWalker = NULL;
3747 double *matrix = ReturnFullMatrixforSymmetric(cell_size);
3748 enum Shading *ColorList = NULL;
3749 double tmp;
3750 vector TranslationVector;
3751 //class StackClass<atom *> *CompStack = NULL;
3752 class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
3753 bool flag = true;
3754
3755 *out << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl;
3756
3757 ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");
3758 while (flag) {
3759 // remove bonds that are beyond bonddistance
3760 for(int i=NDIM;i--;)
3761 TranslationVector.x[i] = 0.;
3762 // scan all bonds
3763 Binder = first;
3764 flag = false;
3765 while ((!flag) && (Binder->next != last)) {
3766 Binder = Binder->next;
3767 for (int i=NDIM;i--;) {
3768 tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
3769 //*out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
3770 if (tmp > BondDistance) {
3771 OtherBinder = Binder->next; // note down binding partner for later re-insertion
3772 unlink(Binder); // unlink bond
3773 //*out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
3774 flag = true;
3775 break;
3776 }
3777 }
3778 }
3779 if (flag) {
3780 // create translation vector from their periodically modified distance
3781 for (int i=NDIM;i--;) {
3782 tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
3783 if (fabs(tmp) > BondDistance)
3784 TranslationVector.x[i] = (tmp < 0) ? +1. : -1.;
3785 }
3786 TranslationVector.MatrixMultiplication(matrix);
3787 *out << Verbose(3) << "Translation vector is ";
3788 TranslationVector.Output(out);
3789 *out << endl;
3790 // apply to all atoms of first component via BFS
3791 for (int i=AtomCount;i--;)
3792 ColorList[i] = white;
3793 AtomStack->Push(Binder->leftatom);
3794 while (!AtomStack->IsEmpty()) {
3795 Walker = AtomStack->PopFirst();
3796 //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
3797 ColorList[Walker->nr] = black; // mark as explored
3798 Walker->x.AddVector(&TranslationVector); // translate
3799 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through all binding partners
3800 if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
3801 OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
3802 if (ColorList[OtherWalker->nr] == white) {
3803 AtomStack->Push(OtherWalker); // push if yet unexplored
3804 }
3805 }
3806 }
3807 }
3808 // re-add bond
3809 link(Binder, OtherBinder);
3810 } else {
3811 *out << Verbose(3) << "No corrections for this fragment." << endl;
3812 }
3813 //delete(CompStack);
3814 }
3815
3816 // free allocated space from ReturnFullMatrixforSymmetric()
3817 delete(AtomStack);
3818 Free((void **)&ColorList, "molecule::ScanForPeriodicCorrection: *ColorList");
3819 Free((void **)&matrix, "molecule::ScanForPeriodicCorrection: *matrix");
3820 *out << Verbose(2) << "End of ScanForPeriodicCorrection." << endl;
3821};
3822
3823/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
3824 * \param *symm 6-dim array of unique symmetric matrix components
3825 * \return allocated NDIM*NDIM array with the symmetric matrix
3826 */
3827double * molecule::ReturnFullMatrixforSymmetric(double *symm)
3828{
3829 double *matrix = (double *) Malloc(sizeof(double)*NDIM*NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
3830 matrix[0] = symm[0];
3831 matrix[1] = symm[1];
3832 matrix[2] = symm[3];
3833 matrix[3] = symm[1];
3834 matrix[4] = symm[2];
3835 matrix[5] = symm[4];
3836 matrix[6] = symm[3];
3837 matrix[7] = symm[4];
3838 matrix[8] = symm[5];
3839 return matrix;
3840};
3841
3842bool KeyCompare::operator() (const KeySet SubgraphA, const KeySet SubgraphB) const
3843{
3844 //cout << "my check is used." << endl;
3845 if (SubgraphA.size() < SubgraphB.size()) {
3846 return true;
3847 } else {
3848 if (SubgraphA.size() > SubgraphB.size()) {
3849 return false;
3850 } else {
3851 KeySet::iterator IteratorA = SubgraphA.begin();
3852 KeySet::iterator IteratorB = SubgraphB.begin();
3853 while ((IteratorA != SubgraphA.end()) && (IteratorB != SubgraphB.end())) {
3854 if ((*IteratorA) < (*IteratorB))
3855 return true;
3856 else if ((*IteratorA) > (*IteratorB)) {
3857 return false;
3858 } // else, go on to next index
3859 IteratorA++;
3860 IteratorB++;
3861 } // end of while loop
3862 }// end of check in case of equal sizes
3863 }
3864 return false; // if we reach this point, they are equal
3865};
3866
3867//bool operator < (KeySet SubgraphA, KeySet SubgraphB)
3868//{
3869// return KeyCompare(SubgraphA, SubgraphB);
3870//};
3871
3872/** Checking whether KeySet is not already present in Graph, if so just adds factor.
3873 * \param *out output stream for debugging
3874 * \param &set KeySet to insert
3875 * \param &graph Graph to insert into
3876 * \param *counter pointer to unique fragment count
3877 * \param factor energy factor for the fragment
3878 */
3879inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment)
3880{
3881 GraphTestPair testGraphInsert;
3882
3883 testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,Fragment->TEFactor))); // store fragment number and current factor
3884 if (testGraphInsert.second) {
3885 *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " successfully inserted." << endl;
3886 Fragment->FragmentCounter++;
3887 } else {
3888 *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
3889 ((*(testGraphInsert.first)).second).second += Fragment->TEFactor; // increase the "created" counter
3890 *out << Verbose(2) << "New factor is " << ((*(testGraphInsert.first)).second).second << "." << endl;
3891 }
3892};
3893//void inline InsertIntoGraph(ofstream *out, KeyStack &stack, Graph &graph, int *counter, double factor)
3894//{
3895// // copy stack contents to set and call overloaded function again
3896// KeySet set;
3897// for(KeyStack::iterator runner = stack.begin(); runner != stack.begin(); runner++)
3898// set.insert((*runner));
3899// InsertIntoGraph(out, set, graph, counter, factor);
3900//};
3901
3902/** Inserts each KeySet in \a graph2 into \a graph1.
3903 * \param *out output stream for debugging
3904 * \param graph1 first (dest) graph
3905 * \param graph2 second (source) graph
3906 * \param *counter keyset counter that gets increased
3907 */
3908inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
3909{
3910 GraphTestPair testGraphInsert;
3911
3912 for(Graph::iterator runner = graph2.begin(); runner != graph2.end(); runner++) {
3913 testGraphInsert = graph1.insert(GraphPair ((*runner).first,pair<int,double>((*counter)++,((*runner).second).second))); // store fragment number and current factor
3914 if (testGraphInsert.second) {
3915 *out << Verbose(2) << "KeySet " << (*counter)-1 << " successfully inserted." << endl;
3916 } else {
3917 *out << Verbose(2) << "KeySet " << (*counter)-1 << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
3918 ((*(testGraphInsert.first)).second).second += (*runner).second.second;
3919 *out << Verbose(2) << "New factor is " << (*(testGraphInsert.first)).second.second << "." << endl;
3920 }
3921 }
3922};
3923
3924
3925/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
3926 * -# constructs a complete keyset of the molecule
3927 * -# In a loop over all possible roots from the given rootstack
3928 * -# increases order of root site
3929 * -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
3930 * -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
3931as the restricted one and each site in the set as the root)
3932 * -# these are merged into a fragment list of keysets
3933 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
3934 * Important only is that we create all fragments, it is not important if we create them more than once
3935 * as these copies are filtered out via use of the hash table (KeySet).
3936 * \param *out output stream for debugging
3937 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
3938 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
3939 * \param *MinimumRingSize minimum ring size for each atom (molecule::Atomcount)
3940 * \return pointer to Graph list
3941 */
3942void molecule::FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
3943{
3944 Graph ***FragmentLowerOrdersList = NULL;
3945 int NumLevels, NumMolecules, TotalNumMolecules = 0, *NumMoleculesOfOrder = NULL;
3946 int counter = 0, Order;
3947 int UpgradeCount = RootStack.size();
3948 KeyStack FragmentRootStack;
3949 int RootKeyNr, RootNr;
3950 struct UniqueFragments FragmentSearch;
3951
3952 *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
3953
3954 // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
3955 // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
3956 NumMoleculesOfOrder = (int *) Malloc(sizeof(int)*UpgradeCount, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
3957 FragmentLowerOrdersList = (Graph ***) Malloc(sizeof(Graph **)*UpgradeCount, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
3958
3959 // initialise the fragments structure
3960 FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
3961 FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
3962 FragmentSearch.FragmentCounter = 0;
3963 FragmentSearch.FragmentSet = new KeySet;
3964 FragmentSearch.Root = FindAtom(RootKeyNr);
3965 for (int i=AtomCount;i--;) {
3966 FragmentSearch.Labels[i] = -1;
3967 FragmentSearch.ShortestPathList[i] = -1;
3968 }
3969
3970 // Construct the complete KeySet which we need for topmost level only (but for all Roots)
3971 atom *Walker = start;
3972 KeySet CompleteMolecule;
3973 while (Walker->next != end) {
3974 Walker = Walker->next;
3975 CompleteMolecule.insert(Walker->GetTrueFather()->nr);
3976 }
3977
3978 // 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
3979 // 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),
3980 // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
3981 // 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)
3982 RootNr = 0; // counts through the roots in RootStack
3983 while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
3984 RootKeyNr = RootStack.front();
3985 RootStack.pop_front();
3986 Walker = FindAtom(RootKeyNr);
3987 // check cyclic lengths
3988 if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
3989 *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;
3990 } else {
3991 // increase adaptive order by one
3992 Walker->GetTrueFather()->AdaptiveOrder++;
3993 Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
3994
3995 // initialise Order-dependent entries of UniqueFragments structure
3996 FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
3997 FragmentSearch.BondsPerSPCount = (int *) Malloc(sizeof(int)*Order, "molecule::PowerSetGenerator: *BondsPerSPCount");
3998 for (int i=Order;i--;) {
3999 FragmentSearch.BondsPerSPList[2*i] = new bond(); // start node
4000 FragmentSearch.BondsPerSPList[2*i+1] = new bond(); // end node
4001 FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1]; // intertwine these two
4002 FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
4003 FragmentSearch.BondsPerSPCount[i] = 0;
4004 }
4005
4006 // allocate memory for all lower level orders in this 1D-array of ptrs
4007 NumLevels = 1 << (Order-1); // (int)pow(2,Order);
4008 FragmentLowerOrdersList[RootNr] = (Graph **) Malloc(sizeof(Graph *)*NumLevels, "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
4009
4010 // create top order where nothing is reduced
4011 *out << Verbose(0) << "==============================================================================================================" << endl;
4012 *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr-1) << " Roots remaining." << endl;
4013
4014 // Create list of Graphs of current Bond Order (i.e. F_{ij})
4015 FragmentLowerOrdersList[RootNr][0] = new Graph;
4016 FragmentSearch.TEFactor = 1.;
4017 FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0]; // set to insertion graph
4018 FragmentSearch.Root = Walker;
4019 NumMoleculesOfOrder[RootNr] = PowerSetGenerator(out, Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
4020 *out << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
4021 NumMolecules = 0;
4022
4023 if ((NumLevels >> 1) > 0) {
4024 // create lower order fragments
4025 *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
4026 Order = Walker->AdaptiveOrder;
4027 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)
4028 // step down to next order at (virtual) boundary of powers of 2 in array
4029 while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
4030 Order--;
4031 *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
4032 for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
4033 int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
4034 *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
4035 *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
4036
4037 // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
4038 //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
4039 //NumMolecules = 0;
4040 FragmentLowerOrdersList[RootNr][dest] = new Graph;
4041 for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
4042 for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
4043 Graph TempFragmentList;
4044 FragmentSearch.TEFactor = -(*runner).second.second;
4045 FragmentSearch.Leaflet = &TempFragmentList; // set to insertion graph
4046 FragmentSearch.Root = FindAtom(*sprinter);
4047 NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
4048 // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
4049 *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
4050 InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
4051 }
4052 }
4053 *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
4054 }
4055 }
4056 }
4057 // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
4058 //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
4059 TotalNumMolecules += NumMoleculesOfOrder[RootNr];
4060 *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
4061 RootStack.push_back(RootKeyNr); // put back on stack
4062 RootNr++;
4063
4064 // free Order-dependent entries of UniqueFragments structure for next loop cycle
4065 Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
4066 for (int i=Order;i--;) {
4067 delete(FragmentSearch.BondsPerSPList[2*i]);
4068 delete(FragmentSearch.BondsPerSPList[2*i+1]);
4069 }
4070 Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
4071 }
4072 }
4073 *out << Verbose(0) << "==============================================================================================================" << endl;
4074 *out << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl;
4075 *out << Verbose(0) << "==============================================================================================================" << endl;
4076
4077 // cleanup FragmentSearch structure
4078 Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
4079 Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
4080 delete(FragmentSearch.FragmentSet);
4081
4082 // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
4083 // 5433222211111111
4084 // 43221111
4085 // 3211
4086 // 21
4087 // 1
4088
4089 // Subsequently, we combine all into a single list (FragmentList)
4090
4091 *out << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl;
4092 if (FragmentList == NULL) {
4093 FragmentList = new Graph;
4094 counter = 0;
4095 } else {
4096 counter = FragmentList->size();
4097 }
4098 RootNr = 0;
4099 while (!RootStack.empty()) {
4100 RootKeyNr = RootStack.front();
4101 RootStack.pop_front();
4102 Walker = FindAtom(RootKeyNr);
4103 NumLevels = 1 << (Walker->AdaptiveOrder - 1);
4104 for(int i=0;i<NumLevels;i++) {
4105 InsertGraphIntoGraph(out, *FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
4106 delete(FragmentLowerOrdersList[RootNr][i]);
4107 }
4108 Free((void **)&FragmentLowerOrdersList[RootNr], "molecule::FragmentBOSSANOVA: **FragmentLowerOrdersList[]");
4109 RootNr++;
4110 }
4111 Free((void **)&FragmentLowerOrdersList, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
4112 Free((void **)&NumMoleculesOfOrder, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
4113
4114 *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
4115};
4116
4117/** Comparison function for GSL heapsort on distances in two molecules.
4118 * \param *a
4119 * \param *b
4120 * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
4121 */
4122inline int CompareDoubles (const void * a, const void * b)
4123{
4124 if (*(double *)a > *(double *)b)
4125 return -1;
4126 else if (*(double *)a < *(double *)b)
4127 return 1;
4128 else
4129 return 0;
4130};
4131
4132/** Determines whether two molecules actually contain the same atoms and coordination.
4133 * \param *out output stream for debugging
4134 * \param *OtherMolecule the molecule to compare this one to
4135 * \param threshold upper limit of difference when comparing the coordination.
4136 * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
4137 */
4138int * molecule::IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold)
4139{
4140 int flag;
4141 double *Distances = NULL, *OtherDistances = NULL;
4142 vector CenterOfGravity, OtherCenterOfGravity;
4143 size_t *PermMap = NULL, *OtherPermMap = NULL;
4144 int *PermutationMap = NULL;
4145 atom *Walker = NULL;
4146 bool result = true; // status of comparison
4147
4148 *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
4149 /// first count both their atoms and elements and update lists thereby ...
4150 //*out << Verbose(0) << "Counting atoms, updating list" << endl;
4151 CountAtoms(out);
4152 OtherMolecule->CountAtoms(out);
4153 CountElements();
4154 OtherMolecule->CountElements();
4155
4156 /// ... and compare:
4157 /// -# AtomCount
4158 if (result) {
4159 if (AtomCount != OtherMolecule->AtomCount) {
4160 *out << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
4161 result = false;
4162 } else *out << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
4163 }
4164 /// -# ElementCount
4165 if (result) {
4166 if (ElementCount != OtherMolecule->ElementCount) {
4167 *out << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
4168 result = false;
4169 } else *out << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
4170 }
4171 /// -# ElementsInMolecule
4172 if (result) {
4173 for (flag=MAX_ELEMENTS;flag--;) {
4174 //*out << Verbose(5) << "Element " << flag << ": " << ElementsInMolecule[flag] << " <-> " << OtherMolecule->ElementsInMolecule[flag] << "." << endl;
4175 if (ElementsInMolecule[flag] != OtherMolecule->ElementsInMolecule[flag])
4176 break;
4177 }
4178 if (flag < MAX_ELEMENTS) {
4179 *out << Verbose(4) << "ElementsInMolecule don't match." << endl;
4180 result = false;
4181 } else *out << Verbose(4) << "ElementsInMolecule match." << endl;
4182 }
4183 /// then determine and compare center of gravity for each molecule ...
4184 if (result) {
4185 *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
4186 DetermineCenter(CenterOfGravity);
4187 OtherMolecule->DetermineCenter(OtherCenterOfGravity);
4188 *out << Verbose(5) << "Center of Gravity: ";
4189 CenterOfGravity.Output(out);
4190 *out << endl << Verbose(5) << "Other Center of Gravity: ";
4191 OtherCenterOfGravity.Output(out);
4192 *out << endl;
4193 if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
4194 *out << Verbose(4) << "Centers of gravity don't match." << endl;
4195 result = false;
4196 }
4197 }
4198
4199 /// ... then make a list with the euclidian distance to this center for each atom of both molecules
4200 if (result) {
4201 *out << Verbose(5) << "Calculating distances" << endl;
4202 Distances = (double *) Malloc(sizeof(double)*AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
4203 OtherDistances = (double *) Malloc(sizeof(double)*AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
4204 Walker = start;
4205 while (Walker->next != end) {
4206 Walker = Walker->next;
4207 Distances[Walker->nr] = CenterOfGravity.Distance(&Walker->x);
4208 }
4209 Walker = OtherMolecule->start;
4210 while (Walker->next != OtherMolecule->end) {
4211 Walker = Walker->next;
4212 OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
4213 }
4214
4215 /// ... sort each list (using heapsort (o(N log N)) from GSL)
4216 *out << Verbose(5) << "Sorting distances" << endl;
4217 PermMap = (size_t *) Malloc(sizeof(size_t)*AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
4218 OtherPermMap = (size_t *) Malloc(sizeof(size_t)*AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
4219 gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
4220 gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
4221 PermutationMap = (int *) Malloc(sizeof(int)*AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
4222 *out << Verbose(5) << "Combining Permutation Maps" << endl;
4223 for(int i=AtomCount;i--;)
4224 PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
4225
4226 /// ... and compare them step by step, whether the difference is individiually(!) below \a threshold for all
4227 *out << Verbose(4) << "Comparing distances" << endl;
4228 flag = 0;
4229 for (int i=0;i<AtomCount;i++) {
4230 *out << Verbose(5) << "Distances: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl;
4231 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold)
4232 flag = 1;
4233 }
4234 Free((void **)&PermMap, "molecule::IsEqualToWithinThreshold: *PermMap");
4235 Free((void **)&OtherPermMap, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
4236
4237 /// free memory
4238 Free((void **)&Distances, "molecule::IsEqualToWithinThreshold: Distances");
4239 Free((void **)&OtherDistances, "molecule::IsEqualToWithinThreshold: OtherDistances");
4240 if (flag) { // if not equal
4241 Free((void **)&PermutationMap, "molecule::IsEqualToWithinThreshold: *PermutationMap");
4242 result = false;
4243 }
4244 }
4245 /// return pointer to map if all distances were below \a threshold
4246 *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
4247 if (result) {
4248 *out << Verbose(3) << "Result: Equal." << endl;
4249 return PermutationMap;
4250 } else {
4251 *out << Verbose(3) << "Result: Not equal." << endl;
4252 return NULL;
4253 }
4254};
4255
4256/** Returns an index map for two father-son-molecules.
4257 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
4258 * \param *out output stream for debugging
4259 * \param *OtherMolecule corresponding molecule with fathers
4260 * \return allocated map of size molecule::AtomCount with map
4261 * \todo make this with a good sort O(n), not O(n^2)
4262 */
4263int * molecule::GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule)
4264{
4265 atom *Walker = NULL, *OtherWalker = NULL;
4266 *out << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
4267 int *AtomicMap = (int *) Malloc(sizeof(int)*AtomCount, "molecule::GetAtomicMap: *AtomicMap"); //Calloc
4268 for (int i=AtomCount;i--;)
4269 AtomicMap[i] = -1;
4270 if (OtherMolecule == this) { // same molecule
4271 for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
4272 AtomicMap[i] = i;
4273 *out << Verbose(4) << "Map is trivial." << endl;
4274 } else {
4275 *out << Verbose(4) << "Map is ";
4276 Walker = start;
4277 while (Walker->next != end) {
4278 Walker = Walker->next;
4279 if (Walker->father == NULL) {
4280 AtomicMap[Walker->nr] = -2;
4281 } else {
4282 OtherWalker = OtherMolecule->start;
4283 while (OtherWalker->next != OtherMolecule->end) {
4284 OtherWalker = OtherWalker->next;
4285 //for (int i=0;i<AtomCount;i++) { // search atom
4286 //for (int j=0;j<OtherMolecule->AtomCount;j++) {
4287 //*out << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
4288 if (Walker->father == OtherWalker)
4289 AtomicMap[Walker->nr] = OtherWalker->nr;
4290 }
4291 }
4292 *out << AtomicMap[Walker->nr] << "\t";
4293 }
4294 *out << endl;
4295 }
4296 *out << Verbose(3) << "End of GetFatherAtomicMap." << endl;
4297 return AtomicMap;
4298};
4299
Note: See TracBrowser for help on using the repository browser.