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