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