source: src/molecules.cpp@ 6ac7ee

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

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

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