source: src/molecules.cpp@ d1df9b

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

new class ConfigFileBuffer and new overloaded function ParseForParaneter() based on this buffer

  • the problem is that we need to use molecuilder for the periodic translating of atoms in a xzy file. However, molecuilder resorts them per element. This destroty the id mapping needed for the DBOND file. The storing was corrected, but also the loading relies on this order.
  • To solve this, we have to pull off quite something: Parse file into buffer, resort the Ion_Type ones (via a map only), load atoms and then put them into the molecule in the original order!
  • Hence, in config::Load we also have a LinearList (map<int, atom*>) in order to add the atoms after the scanning in their original order into the molecule.
  • class ConfigFileBuffer parses a config file line-wise into a buffer and allows for resorting (via a map) of lines containing key_word Ion_Type.
  • BUGFIX: molecule::Output...() were writing the wrong elements, as they were numbered during the output and not before in ascending order! (I.e. if the first atom is Si, but there is also H present, then this Si would become Ion_Type1_1 instead of Ion_Type2_1, because the elements are still sorted by their Z value, hence Ion_Type1 is H!)

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