Changeset a7b761b for src/molecule_fragmentation.cpp
- Timestamp:
- May 27, 2010, 10:46:54 AM (15 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:
- 1024cb
- Parents:
- 8f215d (diff), 05a97c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the(diff)
links above to see all the changes relative to each parent. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/molecule_fragmentation.cpp
r8f215d ra7b761b 39 39 int FragmentCount; 40 40 // get maximum bond degree 41 atom *Walker = start; 42 while (Walker->next != end) { 43 Walker = Walker->next; 44 c = (Walker->ListOfBonds.size() > c) ? Walker->ListOfBonds.size() : c; 41 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 42 c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c; 45 43 } 46 44 FragmentCount = NoNonHydrogen*(1 << (c*order)); … … 359 357 map<double, pair<int,int> > * ReMapAdaptiveCriteriaListToValue(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol) 360 358 { 361 atom *Walker = mol->start;359 atom *Walker = NULL; 362 360 map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ; 363 361 DoLog(1) && (Log() << Verbose(1) << "Root candidate list is: " << endl); … … 391 389 bool MarkUpdateCandidates(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol) 392 390 { 393 atom *Walker = mol->start;391 atom *Walker = NULL; 394 392 int No = -1; 395 393 bool status = false; … … 435 433 bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path) 436 434 { 437 atom *Walker = start;438 435 bool status = false; 439 436 440 437 // initialize mask list 441 for(int i= AtomCount;i--;)438 for(int i=getAtomCount();i--;) 442 439 AtomMask[i] = false; 443 440 444 441 if (Order < 0) { // adaptive increase of BondOrder per site 445 if (AtomMask[ AtomCount] == true) // break after one step442 if (AtomMask[getAtomCount()] == true) // break after one step 446 443 return false; 447 444 … … 457 454 if (AdaptiveCriteriaList->empty()) { 458 455 DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl); 459 while (Walker->next != end) { 460 Walker = Walker->next; 456 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 461 457 #ifdef ADDHYDROGEN 462 if ( Walker->type->Z != 1) // skip hydrogen458 if ((*iter)->type->Z != 1) // skip hydrogen 463 459 #endif 464 460 { 465 AtomMask[ Walker->nr] = true; // include all (non-hydrogen) atoms461 AtomMask[(*iter)->nr] = true; // include all (non-hydrogen) atoms 466 462 status = true; 467 463 } … … 478 474 Free(&FinalRootCandidates); 479 475 } else { // global increase of Bond Order 480 while (Walker->next != end) { 481 Walker = Walker->next; 476 for(molecule::const_iterator iter = begin(); iter != end(); ++iter) { 482 477 #ifdef ADDHYDROGEN 483 if ( Walker->type->Z != 1) // skip hydrogen478 if ((*iter)->type->Z != 1) // skip hydrogen 484 479 #endif 485 480 { 486 AtomMask[ Walker->nr] = true; // include all (non-hydrogen) atoms487 if ((Order != 0) && ( Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))481 AtomMask[(*iter)->nr] = true; // include all (non-hydrogen) atoms 482 if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->nr])) 488 483 status = true; 489 484 } 490 485 } 491 if (( Order == 0) && (AtomMask[AtomCount] == false)) // single stepping, just check486 if ((!Order) && (!AtomMask[getAtomCount()])) // single stepping, just check 492 487 status = true; 493 488 … … 500 495 } 501 496 502 PrintAtomMask(AtomMask, AtomCount); // for debugging497 PrintAtomMask(AtomMask, getAtomCount()); // for debugging 503 498 504 499 return status; … … 516 511 return false; 517 512 } 518 SortIndex = Malloc<int>( AtomCount, "molecule::CreateMappingLabelsToConfigSequence: *SortIndex");519 for(int i= AtomCount;i--;)513 SortIndex = Malloc<int>(getAtomCount(), "molecule::CreateMappingLabelsToConfigSequence: *SortIndex"); 514 for(int i=getAtomCount();i--;) 520 515 SortIndex[i] = -1; 521 516 … … 524 519 525 520 return true; 521 }; 522 523 524 525 /** Creates a lookup table for true father's Atom::Nr -> atom ptr. 526 * \param *start begin of list (STL iterator, i.e. first item) 527 * \paran *end end of list (STL iterator, i.e. one past last item) 528 * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start) 529 * \param count optional predetermined size for table (otherwise we set the count to highest true father id) 530 * \return true - success, false - failure 531 */ 532 bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count) 533 { 534 bool status = true; 535 int AtomNo; 536 537 if (LookupTable != NULL) { 538 Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl; 539 return false; 540 } 541 542 // count them 543 if (count == 0) { 544 for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron 545 count = (count < (*iter)->GetTrueFather()->nr) ? (*iter)->GetTrueFather()->nr : count; 546 } 547 } 548 if (count <= 0) { 549 Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl; 550 return false; 551 } 552 553 // allocate and fill 554 LookupTable = Calloc<atom *>(count, "CreateFatherLookupTable - **LookupTable"); 555 if (LookupTable == NULL) { 556 eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl; 557 performCriticalExit(); 558 status = false; 559 } else { 560 for (molecule::iterator iter = begin(); iter != end(); ++iter) { 561 AtomNo = (*iter)->GetTrueFather()->nr; 562 if ((AtomNo >= 0) && (AtomNo < count)) { 563 //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl; 564 LookupTable[AtomNo] = (*iter); 565 } else { 566 Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl; 567 status = false; 568 break; 569 } 570 } 571 } 572 573 return status; 526 574 }; 527 575 … … 548 596 MoleculeListClass *BondFragments = NULL; 549 597 int *SortIndex = NULL; 550 int *MinimumRingSize = new int[ AtomCount];598 int *MinimumRingSize = new int[getAtomCount()]; 551 599 int FragmentCounter; 552 600 MoleculeLeafClass *MolecularWalker = NULL; … … 576 624 577 625 // create lookup table for Atom::nr 578 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable( start, end, ListOfAtoms, AtomCount);626 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, getAtomCount()); 579 627 580 628 // === compare it with adjacency file === … … 586 634 587 635 // analysis of the cycles (print rings, get minimum cycle length) for each subgraph 588 for(int i= AtomCount;i--;)589 MinimumRingSize[i] = AtomCount;636 for(int i=getAtomCount();i--;) 637 MinimumRingSize[i] = getAtomCount(); 590 638 MolecularWalker = Subgraphs; 591 639 FragmentCounter = 0; … … 624 672 // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle 625 673 KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()]; 626 AtomMask = new bool[ AtomCount+1];627 AtomMask[ AtomCount] = false;674 AtomMask = new bool[getAtomCount()+1]; 675 AtomMask[getAtomCount()] = false; 628 676 FragmentationToDo = false; // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards 629 677 while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) { 630 678 FragmentationToDo = FragmentationToDo || CheckOrder; 631 AtomMask[ AtomCount] = true; // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()679 AtomMask[getAtomCount()] = true; // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite() 632 680 // ===== 6b. fill RootStack for each subgraph (second adaptivity check) ===== 633 681 Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0)); … … 768 816 bool molecule::ParseOrderAtSiteFromFile(char *path) 769 817 { 770 unsigned char *OrderArray = Calloc<unsigned char>( AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");771 bool *MaxArray = Calloc<bool>( AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");818 unsigned char *OrderArray = Calloc<unsigned char>(getAtomCount(), "molecule::ParseOrderAtSiteFromFile - *OrderArray"); 819 bool *MaxArray = Calloc<bool>(getAtomCount(), "molecule::ParseOrderAtSiteFromFile - *MaxArray"); 772 820 bool status; 773 821 int AtomNr, value; … … 872 920 atom *OtherFather = NULL; 873 921 atom *FatherOfRunner = NULL; 874 Leaf->CountAtoms(); 875 876 atom *Runner = Leaf->start; 877 while (Runner->next != Leaf->end) { 878 Runner = Runner->next; 922 923 #ifdef ADDHYDROGEN 924 molecule::const_iterator runner; 925 #endif 926 // we increment the iter just before skipping the hydrogen 927 for (molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end();) { 879 928 LonelyFlag = true; 880 FatherOfRunner = Runner->father; 929 FatherOfRunner = (*iter)->father; 930 ASSERT(FatherOfRunner,"Atom without father found"); 881 931 if (SonList[FatherOfRunner->nr] != NULL) { // check if this, our father, is present in list 882 932 // create all bonds … … 889 939 // Log() << Verbose(3) << "Adding Bond: "; 890 940 // Log() << Verbose(0) << 891 Leaf->AddBond( Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree);941 Leaf->AddBond((*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree); 892 942 // Log() << Verbose(0) << "." << endl; 893 //NumBonds[ Runner->nr]++;943 //NumBonds[(*iter)->nr]++; 894 944 } else { 895 945 // Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl; … … 899 949 // Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl; 900 950 #ifdef ADDHYDROGEN 901 //Log() << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;902 if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), Runner, FatherOfRunner, OtherFather, IsAngstroem))951 //Log() << Verbose(3) << "Adding Hydrogen to " << (*iter)->Name << " and a bond in between." << endl; 952 if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem)) 903 953 exit(1); 904 954 #endif 905 //NumBonds[ Runner->nr] += Binder->BondDegree;955 //NumBonds[(*iter)->nr] += Binder->BondDegree; 906 956 } 907 957 } 908 958 } else { 909 DoeLog(1) && (eLog()<< Verbose(1) << "Son " << Runner->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl); 910 } 911 if ((LonelyFlag) && (Leaf->AtomCount > 1)) { 912 DoLog(0) && (Log() << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl); 913 } 959 DoeLog(1) && (eLog()<< Verbose(1) << "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl); 960 } 961 if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) { 962 DoLog(0) && (Log() << Verbose(0) << **iter << "has got bonds only to hydrogens!" << endl); 963 } 964 ++iter; 914 965 #ifdef ADDHYDROGEN 915 while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen 916 Runner = Runner->next; 966 while ((iter != Leaf->end()) && ((*iter)->type->Z == 1)){ // skip added hydrogen 967 iter++; 968 } 917 969 #endif 918 970 } … … 929 981 molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem) 930 982 { 931 atom **SonList = Calloc<atom*>( AtomCount, "molecule::StoreFragmentFromStack: **SonList");983 atom **SonList = Calloc<atom*>(getAtomCount(), "molecule::StoreFragmentFromStack: **SonList"); 932 984 molecule *Leaf = World::getInstance().createMolecule(); 933 985 … … 1556 1608 FragmentSearch.FragmentSet = new KeySet; 1557 1609 FragmentSearch.Root = FindAtom(RootKeyNr); 1558 FragmentSearch.ShortestPathList = Malloc<int>( AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");1559 for (int i= AtomCount;i--;) {1610 FragmentSearch.ShortestPathList = Malloc<int>(getAtomCount(), "molecule::PowerSetGenerator: *ShortestPathList"); 1611 for (int i=getAtomCount();i--;) { 1560 1612 FragmentSearch.ShortestPathList[i] = -1; 1561 1613 } 1562 1614 1563 1615 // Construct the complete KeySet which we need for topmost level only (but for all Roots) 1564 atom *Walker = start;1565 1616 KeySet CompleteMolecule; 1566 while (Walker->next != end) { 1567 Walker = Walker->next; 1568 CompleteMolecule.insert(Walker->GetTrueFather()->nr); 1617 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 1618 CompleteMolecule.insert((*iter)->GetTrueFather()->nr); 1569 1619 } 1570 1620 … … 1577 1627 RootKeyNr = RootStack.front(); 1578 1628 RootStack.pop_front(); 1579 Walker = FindAtom(RootKeyNr);1629 atom *Walker = FindAtom(RootKeyNr); 1580 1630 // check cyclic lengths 1581 1631 //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) { … … 1664 1714 Vector Translationvector; 1665 1715 //class StackClass<atom *> *CompStack = NULL; 1666 class StackClass<atom *> *AtomStack = new StackClass<atom *>( AtomCount);1716 class StackClass<atom *> *AtomStack = new StackClass<atom *>(getAtomCount()); 1667 1717 bool flag = true; 1668 1718 1669 1719 DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl); 1670 1720 1671 ColorList = Calloc<enum Shading>( AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");1721 ColorList = Calloc<enum Shading>(getAtomCount(), "molecule::ScanForPeriodicCorrection: *ColorList"); 1672 1722 while (flag) { 1673 1723 // remove bonds that are beyond bonddistance … … 1701 1751 Log() << Verbose(0) << Translationvector << endl; 1702 1752 // apply to all atoms of first component via BFS 1703 for (int i= AtomCount;i--;)1753 for (int i=getAtomCount();i--;) 1704 1754 ColorList[i] = white; 1705 1755 AtomStack->Push(Binder->leftatom);
Note:
See TracChangeset
for help on using the changeset viewer.