Chapter 2 How Can One Locate the Global Energy Minimum for Hydrogen-Bonded Clusters? SergeyKazachenkoandAjitJ.Thakkar Abstract Animportantprobleminmanyareasofchemistryandphysicsisfinding theglobalenergyminimumonapotentialenergysurface.Thedifficultystemsfrom the exponential increase in the number of local minima with the size of the sys- tem.Anefficientalgorithmtofindtheglobalminimaofwaterclustersisdescribed andtested.Itworkswellforclusterscontaininguptoabout55watermolecules.A generalization to other hydrogen-bonded clusters is outlined. Applications of this algorithmtowaterclustersandmethanolclustershavealreadybeenreportedinthe literature. 2.1 Introduction Material particles consisting of a few to a few thousand atoms are called clusters. Clusterpropertiescanhavedramaticsizeandshapedependence.Clusterscanserve asbuildingblocksfornewmaterialsandelectronicdevices.Henceclustersofmet- als, semiconductors, ionic solids, rare gases, and small molecules have been stud- ied using both theoretical and experimental methods. Atomic and molecular clus- ters[1]areheldtogetherbyhydrogenbonds[2]orbyrelativelyweakintermolecular forces[3, 4]. Inparticular,smallclustersofhydrogen-bondedwatermoleculeshavereceiveda lotofattention;see,forexample,spectroscopicwork[5–16],anddensityfunctional theoryandabinitioinvestigations[17–31].Waterisnotasimplesubstanceandhas anomalousphysicalandchemicalproperties.Morethanacenturyofworkhasbeen devotedtomodelingandunderstandingtheseproperties.Nevertheless,manyaspects of water remain unsolved puzzles and the development of water models continues S.Kazachenko⋅A.J.Thakkar(✉) DepartmentofChemistry,UniversityofNewBrunswick, Fredericton,NBE3B5A3,Canada e-mail:[email protected] S.Kazachenko PresentAddress:DepartmentofChemistry,Queen’sUniversity, Kingston,ONK7L3N6,Canada ©SpringerScience+BusinessMediaNewYork2016 25 J.LeszczynskiandM.K.Shukla(eds.),PracticalAspects ofComputationalChemistryIV,DOI10.1007/978-1-4899-7699-4_2 26 S.KazachenkoandA.J.Thakkar to attract much interest. Water models can use water clusters, (H O) , as building 2 n blocks[32–36]. An important part of cluster research is a characterization of the global min- imum of the potential energy surface—the structure of the most stable cluster. Ascholarlyandcomprehensivesurveyofthedifficultproblemoflocatingtheglobal minimumonapotentialenergysurfaceisavailable[37].Computationalapproaches tothedeterminationofglobalminimaareveryattractivebecausetheexperimental determinationofthestructureofthegroundstateisextremelydifficultforcomplex systems.Unfortunately,thecomputationalsolutionisratherexpensivebecausethe numberoflocalminimaonthesurfacegrowsexponentiallywiththedimensional- ity of the surface or, in other words, with the number of atoms in the system. An algorithm searching for a global minimum can easily get trapped in one of these local minima because the traversal of many intermediate minima and the crossing ofhighbarriersseparatingthemmayberequiredtofindalowerlocalminimum.In high-dimensionalcases,whicharethenormratherthantheexceptioninphysically important clusters, a global minimum can only be discovered and verified after a sufficientlylargenumberoflowenergyminimahavebeenlocatedandcompared. We have been studying water clusters for some time now. Hence we have been gradually developing and refining our own algorithm for locating a global energy minimum for water clusters. We previously described successive improvements to our algorithm in an incremental manner in a series of papers [38–41] and a the- sis [42]. The purpose of this chapter is to try and describe the current state of our algorithminarelativelyself-containedandcohesivemannerwithoutthereaderhav- ingtoconsultourpreviouswork. Allalgorithmtestsonwaterclustersdiscussedinthischapteraredonewiththe empirical TIP4P model [43], a reparametrization of the venerable Bernal-Fowler model [44]. The interactions are considered to be pair-wise additive and the water monomers are held rigid so that they do not vibrate. The interaction between a pairofwatermoleculesisgivenbyaLennard-Jones(12,6)interactionbetweenthe oxygen atoms and electrostatic interactions between three point charges on each water molecule. There are positive charges on the hydrogen atoms and a balanc- ingnegativechargebetweenthehydrogenatomsalongtheC symmetryaxis.The 2 TIP4Pmodelgivesareasonablethermodynamicdescriptionofliquidwater.Many studies have been devoted to finding the global minima of TIP4P water clusters [38–40, 45–51]. Generally good agreement has been found [48] between its pre- dictionsforsmallclusterswithn≤12andbothabinitioandexperimentalresults. TheTIP4Pglobalminimaforclustersizesupto47arenowfirmlyestablished[41] exceptforn=39andn=45.PutativeTIP4Pglobalminimahavebeenreported[41] forlargerclusterswithn≤55.Low-energyTIP4Pstructuresformuchlargerclusters withselectedsizeshavebeenlocated[7, 49].PurewaterclustersbasedontheTIP4P modelarenowabenchmarkformethodsforglobaloptimizationofhydrogen-bonded clusters. Some of the basic terminology commonly used in the global optimization lit- erature is summarized in Sect.2.2 which also contains a brief description of basin andminimahopping.Next,thelongSect.2.3detailsthevariouswaysbywhichwe 2 HowCanOneLocatetheGlobalEnergy... 27 optimize the topology of water clusters with a fixed oxygen skeleton. The overall algorithm is then described in Sect.2.4 and its generalization to other hydrogen- bonded clusters is outlined in Sect.2.5. A few concluding remarks are made in Sect.2.6. 2.2 BasinHoppingandMinimaHopping Adiscussionofsearchalgorithmsisfacilitatedbytheintroductionofsometermi- nologythatiswidelyusedintheliterature.Abasinisaregionofgeometricalconfig- urationspacearoundaminimumonthepotentialenergysurface.Thebasincontains allstructures,orconfigurations,fromwhichasearchcanreachthisminimumusing onlysmallstepsanddownhillmoves.Asuper-basinistheunionofseveralneigh- boringbasins.Afunnel isasuper-basinwiththepropertythatstartingatanypoint in onecanreachthelowestminimumin withoutcrossingbarriersthatarevery highrelativetotheaverageenergydifferencebetweenlocalminimain.Figure2.1 isaschematicillustrationoftwosuper-basins,oneofwhichisafunnelandtheother oneisnot. UseoftheBoltzmannfactor,exp(−ΔE∕kT),tocontrolallthestepsusedtoleave a basin makes the crossing of high barriers a rare event. Hence, global minimiza- tionmethodsbasedexclusivelyonthermodynamicprinciplescanbeextremelyslow becausetheymayfindithardtoexitafunnelandmayrepeatedlyvisitneighboring configurationsthatarecloseinenergy.Manystandardalgorithmssuchassimulated annealing[53–55]andbasinhopping[46,56–59]arebasedonthermodynamicprin- ciples.However,geneticalgorithms[60–62]andminimahopping[52]arenot. The fundamental idea in basin hopping is that the potential energy surface is effectively transformed to a stepped surface [46, 56, 63]. This is done by ending eachsearchstepwithalocaloptimizationsothatoneeffectivelysearchesoverlocal minima.Howonestartsthenextsearchstepdistinguishesvariousalgorithmsbased on this seminal idea. We chose minima hopping [52] as the global optimization Fig.2.1 Twosuper-basins:theoneontherightisafunnelbuttheoneontheleftisnot.Figure adaptedfromGoedecker[52] 28 S.KazachenkoandA.J.Thakkar algorithminourpreliminarywork.Minimahoppingcanbethoughtofasaversion ofbasin-hopping[46, 56, 63].Itsgoalistoavoidrevisitingpreviouslylocatedlocal minimawithoutforbiddingrepeatedpassagesthroughanytransitionbasinsthatmay separatemanyfunnels.Thisisachievedintheminima-hoppingmethodbymaintain- ingalistofvisitedminimaandusingadaptivethresholdstoleaveabasinmoreand morevigorouslyeachtimethealgorithmrevisitsit.Moleculardynamics(MD)steps areusedasanefficientmechanismforcrossingenergybarriers. Theminima-hoppingalgorithmhasaninnerpartforjumpingtoanewlocalmin- imum and an outer part for accepting or rejecting the new local minimum. In the innerpart,onetriestoescapethecurrentminimum𝐌 byashortMDsimulationin c whichtheatomshaveaBoltzmannvelocitydistributionsuchthattheirkineticenergy isfixedatE .Thesimulationisstoppedassoonasaminimumisencounteredalong kin thetrajectoryorthemaximumnumberofMDstepsisexceeded.Thenoneoptimizes totheclosestlocalminimum𝐌usingasuitablelocaloptimizationmethod.Ifthis minimumhasbeenvisitedpreviously,thenmultiplyE by𝛽 >1tomakethenext kin escapeattemptmorevigorousandrepeattheinnerpart.Ifthisisanewlyfoundlocal minimum, divide E by 𝛽 to make the next escape attempt less vigorous and go kin totheouterpart.ThelatteracceptsorrejectsthelocalminimumE(𝐌)asfollows. IfE(𝐌)−E(𝐌 )<E ,theminimumisacceptedandE isdividedby𝛼 >1to c diff diff makethenextacceptancemoredifficult.OtherwiseitisrejectedandE ismulti- diff pliedby𝛼tomakethenextacceptanceeasier.Wetriedvariousvaluesof𝛼and𝛽but wereunabletoimproveuponthevaluesof𝛼 =1.02and𝛽 =1.05recommendedby Goedecker[52].Theminima-hoppingalgorithmstopswhenE reachesorexceeds kin amaximumvalueEmaxorthenumberofinnerstepsexceedsapresetlimitN . kin iter It is important to do local optimizations as efficiently as possible because they areatime-consumingpartofminimahopping.Inourimplementation,weperform local minimization in two steps. The first optimization uses the limited memory L-BFGSmethod[64, 65]withalooseconvergencethreshold.Theoptimizationis refinedinthesecondstepwhichusesDavidon’soptimallyconditionedvariablemet- ricmethod[66]andamorestringentconvergencecriterion. It soon became apparent to us that minima hopping did not always succeed in findingtheoptimumhydrogenbondtopologyforwaterclusters[39].Nevertheless, minima hopping served us as the basic algorithm upon which improved and spe- cializedmethodsforwaterclusterswerebuilt.Specialtopologyrefiningalgorithms, describedinthenextsection,areusedbothwithintheglobalsearchalgorithmand independentlyforarefiningsteponthelistofsavedlowenergyminima.Although manyfeaturesdescribedhereweredesignedtoworkonlywithwater,thecorealgo- rithmisgeneralandhasbeenappliedsuccessfullytootherhydrogen-bondedclusters includingpuremethanolclusters[67]andpureclustersofethanol,n-propanol,and iso-propanol. 2 HowCanOneLocatetheGlobalEnergy... 29 2.3 FindingtheOptimalHydrogenBondTopology 2.3.1 RepresentationofH-BondedClusters Analgorithmforglobaloptimizationofhydrogen-bondedclustersrequiresaconve- nientdescriptionofaclusteranditsfeatures.Atthemostbasiclevel,theclusteris defined by the Cartesian coordinates and types of the atoms it consists of. We use thatinformationtoderivesomeimportantproperties andfindaconvenient wayto representthem. AbasicpropertyofanyH-bondedclusteristhepresenceorabsenceofahydrogen bondbetweenapairofgivenmolecules.Weusesimplegeometriccriteriatodecide whetherthereisanO–H⋯Ohydrogenbond:theH⋯Odistanceshouldbelessthan 2.5ÅandtheO–H⋯Oangleshouldbegreaterthan90◦.Itisconvenienttodescribe theH-bondsinaclusterofnmoleculesbyagraph.Thelatterisrepresentedasan n×nadjacencymatrix𝐀withelementsobeyingtheusualrules:A =1ifthereisa ij bondbetweenmoleculesiandjandA =0otherwise.Bydefinition,A =0because ij ii moleculesarenotconnectedtothemselvesandA =A becausethebonddirection ij ji is not taken into consideration. The adjacency matrix allows one to ca∑lculate, for example,thetotalnumberofhydrogenbondsinaclusterasN = 1 A .One bond 2 ij ij canalsocalculatethenumberofringsofagivensizeformedbyconnectedmolecules. Weimplementedthecountingofringsusingabacktrackingalgorithmbasedonthe workofFranzblau[68].Thereareseveralpossibledefinitionsofaringinagraph.We countallringsinwhicheachmonomerisconnectedtoexactlytwoothermonomers belongingtothesamering. TheH-bonddirectionalityisimportantintopologyoptimization.Keepingtrack ofboththeexistenceanddirectionoftheH-bondsinaclustercanbeaccomplished byusingadigraph.Thelattercanberepresentedbyadirectedadjacencymatrix𝐃. Aswiththeadjacencymatrix,theabsenceofabondbetweenmonomersiandjis indicatedbyD =0andsoalldiagonalelementsvanish,D =0.Ifthereisabond ij ii inthei→jdirectionthenD =1andD =0.The𝐃matrixcontains,forexample, ij ji informationabout∑thenumberof∑donor(don)andacceptor(acc)H-bondsforeach molecule:Ni = D ,Ni = D .Notethatthesummationindexdependson don j ji acc j ij thedefinitionoftheH-bonddirection.Inacomputerprogram,thesparsityof𝐀and 𝐃canbeexploitedbystoringthemaslinkedlistsofnon-zeroelements. 2.3.2 WhyIsTopologyOptimization Needed? Asimplecharacterizationofaclusterofnwatermoleculesisitsskeletonorgraphby whichwemeantheconnectivityofthemonomersdescribedbyanadjacencymatrix 𝐀. The positions of the oxygen atoms and the skeleton define the shape of a clus- ter as in the left panel of Fig.2.2. However, virtually all hydrogen bonds between water molecules can be assigned a direction, say from donor to acceptor. Hence, 30 S.KazachenkoandA.J.Thakkar Fig. 2.2 An example of a framework and a complete cluster with one of the possible H-bond topologiesfora(H O) cluster.Twointernalmoleculesareshowninblue 2 26 there can be many water cluster geometries with the same skeleton but different hydrogen bond topologies; that is to say, cluster structures with the same skeleton butdifferentdirectionsinoneormoreofthehydrogenbonds.Intheterminologyof Sect.2.3.1,anygivenadjacencymatrix𝐀cancorrespondtomanydirectedadjacency matrices𝐃 .TheoxygenframeworkandtheH-bondtopologytogetherdefineacom- k pleteclusterasintherightpanelofFig.2.2.Thestabilityofagivenframeworkcan varysignificantlydependingonthehydrogenbonddistribution[69–75]. Since each cluster framework maybe paired with a large number of H-bond topologies, the resulting cluster structures can have significant energy differences. Therefore,itisimportantthatweareabletolocatetheonewiththelowestenergy.It mightalsobeofinteresttoknowhowmanytopologiesexistforagivenframework andpossiblyseparatethemintocategories.Locatingtheminimumenergytopology, or topology optimization, can be done either as a separate procedure on selected waterclustersorasapartofaglobaloptimization. The problem of finding the best hydrogen bond topology for water clusters has been studied extensively; see, for example, Refs. [49, 51, 69–74, 76–83]. Polyhe- dral and cubic water clusters were described in several publications using graph theory[69, 74, 78, 79].TheeffectsofH-bondtopologyonthestabilityandspec- troscopic properties of water octamers were studied by Francisco et al. [73, 80]. Therelationbetweenthetopologyandinteractionenergywasstudiedforpolyhedral clusters[72, 81]andforsomeothershapesandsizes[71, 81].Asaresult,several formulasweredeveloped topredicttherelativeenergy ofawaterclusterbasedon itsmonomerconnectivity.Adifferentapproachistouseaprotontransfertochange thedirectionofhydrogenbondsandsosampleanumberoftopologiesusinggeneral optimization methods [49, 51]. There is also the brute-force approach of examin- ingallpossibletopologiesforanarbitraryclustershape[77, 82, 83].Thelasttwo approacheswillbediscussedlaterinmoredetail. Unfortunately,mostofthesuggestedmethodsweredesignedonlyforaparticu- laroxygenframeworkorturnedouttobeinefficient.Wecreatedseveralalgorithms suitabletoourgoalofreliablyfindingthelowestenergytopologyforawaterclus- ter of an arbitrary shape. The algorithms were efficient at locating lowest-energy topologieswhenappliedtoclusterswithnomorethan55molecules. 2 HowCanOneLocatetheGlobalEnergy... 31 2.3.3 ComparingandStoringWaterClusterMinima In our approach for H-bond topology optimization, a water cluster is defined as a structure with a particular oxygen framework. A framework is assumed to have a uniqueadjacencymatrix𝐀.Notwoframeworkswillhavethesamematrixand,there- fore,clusterscanbecomparedbycomparingtheir𝐀matrices.Allpossibletopolo- giesforaframeworkareconsideredtobevariationsofthesamecluster.Therefore, eachclusterdefinedbyits𝐀matrixhasasetoflocalminimacorrespondingtoall possibleH-bondtopologies.Suchseparationoftheframeworkandtopologyhelpsto reducethenumberoflocalminimathatarestoredduringaglobalminimumsearch. Onlytheversionofaclusterwiththebesttopologyiskept. Howcanwedetermineiftwosetswithanequalnumberofwatermoleculesform clusterswithanidenticalshape?𝐀matricescanbecompareddirectlyiftheposition of the molecules has changed only slightly, i.e. after a small distortion. However, ingeneralwewantthecomparisontobeindependentoftheorderofmoleculesin a coordinate list. In our case this is done in three steps which helps to reduce the number of time-consuming operations. First, the number of H-bonds must be the same, which is easy to calculate and to check. Next, we use the idea that the way moleculesareconnectedhasaneffectonthenumberofringsformedbyH-bonded monomersinacluster.Moreover,thenumberofringsisindependentoftheorderof molecules in a list.Therefore, as the second step, it is required that the number of ringsofeachsizefrom3to10moleculesmustbethesame.Forthosecaseswhere theringruleisalsosatisfied,analignmentof𝐀matricesisusedtoperformthefinal checkwhichmustallowforthepossibilitythattheorderoftheelements(monomers) is different. A backtracking procedure is used to examine the connectivity of each monomerandtotryandmatchpairsofmoleculesfromtheclustersbeingcompared. Thisprocedureallowsustodeterminewhethertwomatricescorrespondtothesame frameworkevenwhentheorderoftheelementsisnotthesame. Whenthelistofvisitedminimaislargeenough,itisnolongerfeasibletokeep eithertheCartesiancoordinatesortheadjacencymatrixforeachlocalminimum.In thatcase,thethirdstepthatincludesalignmentof𝐀matricescanbeomitted.Nine integers,thenumberofH-bondsandthenumbersofthree-throughten-membered rings, are saved for each local minimum and provide a robust and reliable way to compareclusterframeworks. 2.3.4 FilterstoScreenTopologies Findingthelowest-energytopologyofalargeclusterframeworkwouldbeimpossi- bleifitrequiredthelocaloptimizationofallorevenmostofitstopologies.Fortu- nately,thenumberofexpensivelocaloptimizationsrequiredcanbereducedbyusing simplefiltersorcriteriabasedsolelyongeometricalconsiderationstoweedoutthe topologies that are likely to have a high energy. Note that although the shape of a 32 S.KazachenkoandA.J.Thakkar Fig.2.3 Examplesoflessstableandmorestablehydrogenbondingofamonomer water cluster depends somewhat on its topology, a stable (i.e. low-energy) frame- workdoesnotbreakorchangewhenitstopologyischanged. ThefirstthingtodoistocountthenumberofH-bondsandtheirtypesforeach molecule.Wedefineadonorbondasabondcreatedbyahydrogenatomofagiven moleculeandanacceptorbondastheonecreatedbyanoxygenatomofthemolecule. WeuseamodificationofthevenerableBernal-Fowler“icerules”[44]andrequire that 1≤number of donor bonds≤2, and 1≤number of acceptor bonds≤2or 3. In words, each water molecule in a cluster must have at least one donor and one acceptorbond.Moreover,therecanbeatmosttwodonorbonds(oneforeachhydro- gen)andatmosttwoacceptorbonds.Theserulesallowonetoavoidunstablecon- nectivityofmolecules;seeFig.2.3.Inrarecasesawatermoleculeinaclusteraccepts threehydrogenbondsanditisnecessarytoaccountforsuchapossibility.However, allowingallmoleculestoacceptthreeH-bondswouldsignificantlyincreasethenum- beroftopologycombinations.Thus,thealgorithmdetectspenta-coordinatedmole- culesintheinputgeometryandonlythosecanhaveuptothreeacceptorbonds,while therestarerestrictedtohavenomorethantwoacceptorbonds. We found that in some cases changing the cluster topology leads to unrealistic monomer angles; see Fig.2.4. Note that the molecules in a cluster shift somewhat during a local optimization of a new topology. Thus we must allow for a range of angles to be accepted. Test calculations suggest that a topology can reliably be Fig.2.4 Anexampleofan unrealisticmonomeranglein oneoftheH-bondtopologies foran(H O) cluster 2 12 2 HowCanOneLocatetheGlobalEnergy... 33 Fig. 2.5 Adjacent monomers with a connectivity that leads to higher energy (D—donor, A— acceptor) considered unfavorable if it leads to a monomer angle larger than 150◦ or smaller than65◦.Thesevaluesspanalargeenoughrangetoensurethatnoimportanttopolo- giesaremissedbecauseoftheanglefilter. Itiswell-establishedthatadjacentdanglinghydrogenatomsleadtoahigherclus- terenergy [69, 71–73, 79].Inaddition,Anickfoundtwootherpatterns leadingto anincreasedenergy[81].ThesethreemotifsofconnectivityareshowninFig.2.5. ThetypesofH-bondsarelabeledbyDfordonorandAforacceptor.Thenumberof motifs1and2thatoccurinaclustermustbekeptaslowaspossible.Thismeansthat onecansafelydiscardalltopologieswithmoresuchmotifsthanthelowestnumber foundatanygivenstage.Thethirdruleisnotasstrong.Wefoundthatatopologycan bediscardedsafelyonlywhenthenumberofsuchmotifsexceedsthelowestnumber foundplustwo. 2.3.5 TopologyOptimization byEnumeration A straightforward way of finding the best topology is to generate all possible H- bond distributions, do local optimizations on each one, and compare their ener- gies.MiyakeandAidausedadjacencyanddirectedadjacencymatrices(graphsand digraphs)todescribetheframeworkandthetopologyofacluster,respectively[82]. Knowingthedigraphandthecoordinatesoftheoxygenatoms,onecanrecreatethe topologyandthenperformalocaloptimizationtorelaxtheclustergeometry.Theuse ofadjacencymatriceswasapromisingidea;however,intheiralgorithmallpossible 𝐃matricesweregeneratedfirstandonlylaterweretheycheckedtoseeiftheywere useful. This led to a large wasted computational effort and the largest cluster size they were able to study was limited to just eight water molecules. Vukičević et al. improvedtheperformanceofthemethodbyeliminatingtheneedtogenerateunre- alisticmatrices[77, 83]andanalyzedthetopologyforclusterswithupto12water molecules. Theideasofusinggraphsanddigraphsandeliminatingundesiredcombinations on the fly were used in the creation of our method, called topology enumeration (NT).ConsiderthecagehexamershowninFig.2.6aasanexample.Westartfrom thegivenframeworkanditsadjacencymatrix𝐀.Eachrowandcolumnof𝐀withthe sameindexcorrespondstoamolecule.Definingallconnectionsofamoleculetothe 34 S.KazachenkoandA.J.Thakkar Fig.2.6 Adjacencymatricesfortheenumerationmethod.asamplegeometry;badjacencymatrix; cpossibledirectedadjacencymatrices othermoleculesautomaticallydefinestheconnectionsofallmoleculestothecurrent one. Therefore, it is sufficient to consider each row and column starting from the correspondingdiagonalelement;seeFig.2.6b.Theconnectivityofthelastmolecule (number6),asexpected,iscompletelydefinedbytheconnectivityoftheprevious molecules. Thepositionswithonesineachrow(orcolumn)of𝐀telluswhichmoleculesare connected to a given molecule. Let us choose rows as a reference. Then, by going through all combinations of 1 and 0 for positions marked with “1∕0” in Fig.2.6c wecangenerateallpossibleH-bonddirectionsforthegivenskeleton.Thiscanbe achieved by taking a bit representation of an integer value that ranges from 0 to 2n−1.Fromthedefinitionof𝐃,thevaluesinacolumnaretheoppositeofthevalues inarow(positionswith“0∕1”inFig.2.6c).H-bonddirectionsaregeneratedforeach moleculeinthismannerusingabacktrackingloop.Thebacktrackingallowsoneto applygeometryfiltersforeachmoleculerightawaywithoutgeneratingacomplete matrix.ThestepsconstitutingourbacktrackingmethodareshowninListing2.1. Listing2.1 Backtrackingalgorithmfortopologyenumeration. --------------------------------------------------------------- Mark all bit combinations of all molecules as not used Select molecule 1 as the current molecule ’M’ WHILE (’M’ > 0) IF (all bit combinations for ’M’ have been used) THEN Mark all bit combinations for ’M’ as not used Set ’M’ = ’M’ - 1 ELSE Choose next bit combination for ’M’ Mark this bit combination as used IF (topology filters are satisfied) THEN Set ’M’ = ’M’ + 1 END IF IF (’M’ > number of molecules) THEN Create and save a cluster geometry for local optimization Set ’M’ = ’M’ - 1 END IF END IF END WHILE ---------------------------------------------------------------
Description: