source: src/molecules.cpp@ bf46da

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since bf46da was bf46da, checked in by Frederik Heber <heber@…>, 17 years ago

new function molecule::CreateMappingLabelsToConfigSequence() creates SortIndex in molecule::FragmentMolecule()

further and now final abstraction of FragmentMolecule()

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