source: src/molecules.cpp@ d3a46d

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 d3a46d was d3a46d, checked in by Frederik Heber <heber@…>, 17 years ago

added --enable-optimization, thrown -g -O specifics out of Makefile.am's

--enable-optimization is added to each configure.ac and is linked with --enable-debug (no debug and optimization at the same time, debug is normally off and -O2 on), implemented via C(XX)FLAGS
Henceforth, -g, -O and -cc=... specifics are thrown out of the Makefile.am's, so that the optimiziation and debugging is streamlined throughout the whole package

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