UFDC Home  myUFDC Home  Help 
SOFFT:
EPQCLE: s P SOFFT: s P Figure 4{12. Exp ectation of momen tum and disp ersion for the Nasurface system.
PAGE 72
58 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 5 10 15 20 25 30 35 40 r (R)R (au) EPQCLE SOFFT Figure 4{13. Densit y function ( R ) for the Nasurface system.
PAGE 73
59 system, the coherence rapidly diminishes once the in teraction p oten tial is crossed. This can b e explained b y the size of the energy gap at asymptotic distances, as larger ionization energies lead to more rapidly v anishing coherences. Figure 4{11 sho ws the exp ectation of the p osition of Na steadily increases from its initial a v erage, while the disp ersion in p osition initially decreases and then increases again. On the other hand, Figure 4{12 sho ws a mark ed dierence in the b eha vior of the exp ectation of the momen tum, where it b egins at h P i = 0 au, rapidly increases and then b ecomes stationary around h P i = 78 au. W e also nd that the momen tum disp ersion rst increases, and then decreases. Finally Figure 4{13 presen ts the densit y function at the end of the sim ulation. Because of the distortion of the phase space grid, the densit y function obtained from the EPQCLE had to b e calculated using bins and through appro ximations within these bins, and th us is sub ject to some noise. Nev ertheless, w e nd excellen t agreemen t b et w een the EPQCLE and SOFFT results. F or all observ ables, the EPQCLE results are quan titativ ely similar to the SOFFT results to visual resolution. Ha ving studied b oth Na and Li approac hing a metal surface, w e see that the EPQCLE can b e exp ected to yield v ery accurate results for these kinds of systems, ev en when coherence is main tained o v er long p erio ds. 4.4 Binary Collision In v olving Tw o Av oided Crossings 4.4.1 Mo del Details F or the second system, w e lo ok at a t w ostate collision mo del where the diabatic surfaces in tersect t wice. Because of the dual crossing and the coupling in this region, quan tum in terference and eects suc h as tunnelling pla y a substan tial role in the dynamics of the quasiclassical v ariable. As suc h, this mo del is quite demanding for mixed quan tumclassical metho ds, where one can exp ect deviations at lo w er energies.
PAGE 74
60 T able 4{2. P arameters used in the dual a v oided crossing collision mo del. P arameter V alue (au) R 0 8 P 0 [10, 30] 2.5176 M 2000 The Hamiltonian elemen ts are 79 H 11 ( R ) = 0 ; (4.22) H 22 ( R ) = 0 : 1 exp ( 0 : 28 R 2 ) + 0 : 05 ; (4.23) H 12 ( R ) = 0 : 015 exp( 0 : 06 R 2 ) : (4.24) In the ab o v e, R is the n uclearn uclear separation. The initial state is a ground state Gaussian w a v epac k et whic h b egins in the asymptotic region R 0 = 8 au. Its PWT o v er R giv es the Gaussian densit y 11 ( P ; R ) = 1 R R 0 2 2 ( P P 0 ) 2 # ; (4.25) with 12 = 21 = 22 = 0. A t t = 0, the w a v epac k et propagates to w ard the region of coupling, where it is partially transmitted and partially rerected, no w with p opulations in b oth the ground and excited state. The parameters used in the calculation are sho wn in T able 4{2 The diabatic p oten tials used in the mo del are sho wn in Figure 4{14 4.4.2 Prop erties of In terest As in the alk ali atomsurface mo del, w e can follo w the p opulations o v er time. The system b egins in the ground state, and as it passes through the collision region, w e exp ect some of the p opulation to transfer to the excited state. The amoun t of transfer dep ends on the collision energy W e can also follo w the coherence once the
PAGE 75
61 0.06 0.04 0.02 0 0.02 0.04 0.06 0.08 10 5 0 5 10 Energy (au)R (au) H 11 H 22 H 12 Figure 4{14. P oten tial curv es for the dual a v oided crossing collision.
PAGE 76
62 collision has ended and the system reac hes asymptotic v alues. W e also consider the transmission of the ground state as the w a v epac k et passes through the in teraction. T ransmission. As the w a v epac k et passes through the collision, part of it con tin ues forw ard while the remainder is rerected bac k. W e can study this b y computing the probabilit y of transmission of the ground state, T 1 = 1 Z 0 dR dP 11 ( R ; P ) : (4.26) Although not included, w e could also compute the probabilit y of rerection of the ground state, as w ell as the asymptotic p opulations of the excited state. 4.4.3 Results In Figure 4{15 w e see the total p opulation transfer at energy P 0 = 30 au. In this energy region, the EPQCLE matc hes the quan tum results to within a couple p ercen t. F or b oth the EPQCLE and SOFFT sim ulations, the rerection of the w a v efunction (or PWTDM) w as negligible, so that w e need only compare the transmission. This negligible rerection w as found throughout the energies studied. In Figure 4{16 w e see the coherence as a function of time. The v ariations are large initially and diminish asymptotically to regular oscillations. The EPQCLE captures this b eha vior qualitativ ely and to within 5% quan titativ ely W e sho w the deformation of the phase space grid in Figure 4{17 W e see that it is primarily the inner part of the grid that undergo es deformation, while the enclosing p oin ts main tain a structured order. This is lik ely due to only a small range of v alues in phase space whic h are seriously aected b y the p oten tials. Outside the in teraction region, all HellmannF eynman forces are zero, so one w ould exp ect deformations only for p oin ts whose time in the in teraction region w as signican t. In Figure 4{18 w e displa y the transmission probabilit y for a wide range of energies. As exp ected, the EPQCLE p erforms b etter for higher energies, and while the doublew ell is qualitativ ely repro duced, the EPQCLE fails to giv e quan titativ ely
PAGE 77
63 accurate results for energies lo w er than (log E = 2). This is lik ely due to the quantum tunnelling eects describ ed earlier, whic h are not exp ected to b e repro duced w ell b y the EPQCLE. W e also compare transmission probabilities obtained through surface hopping metho ds b y T ully and co w ork ers. 51 These surface hopping calculations sho w deviations from the quan tal results that are similar to the EPQCLE probabilities at lo w er energies; at higher energies, the EPQCLE is sligh tly sup erior to the surface hopping sc heme. 4.5 Photoinduced Disso ciation of a Diatomic System 4.5.1 Mo del Details F or the third test system, w e explore a mo del of the NaI complex. As in the previous mo dels, it in v olv es t w o diabatic surfaces and an in teraction around the a v oided crossing. A substan tially dieren t feature, ho w ev er, is a longrange Coulom bic attraction in an ionic state. As w e shall see, this attractiv e p oten tial results in the the complex oscillating b et w een a neutral and ionic state as the so dium and io dine separate and come bac k together, partially disso ciating at eac h crossing in to an asymptotically neutral state. This oscillatory motion is a go o d test for the EPQCLE at long times in cases where asymptotic states are not reac hed quic kly The Hamiltonian elemen ts are 86 H 11 ( R ) = A 1 exp [ 1 ( R R 0 )] ; (4.27) H 22 ( R ) = [ A 2 + ( B 2 =R ) 8 ] exp ( R = ) 1 =R ( + + ) = 2 R 4 C 2 =R 6 2 + =R 7 + E 0 ; (4.28) H 12 ( R ) = A 12 exp [ 12 ( R R x ) 2 ] : (4.29) T o form the initial state, w e T a ylor expand (to rst order) the ionic p oten tial H 22 ab out its minim um R = R 0 and nd the lo w est energy state of this harmonic w ell. A t t = 0, the w a v epac k et undergo es a sudden optical promotion to the neutral curv e,
PAGE 78
64 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 200 400 600 800 1000 1200 1400 Populations h1 and h2time (au) EPQCLE: h 1 SOFFT: h 1 EPQCLE: h 2 SOFFT: h 2 Figure 4{15. P opulations 1 and 2 vs. time for the dual crossing collision mo del.
PAGE 79
65 0.2 0.15 0.1 0.05 0 0.05 0.1 0 200 400 600 800 1000 1200 1400 Coherence: Re( h12)time (au) EPQCLE SOFFT Figure 4{16. Coherence describ ed b y R e ( 12 ) vs. time, for the dual crossing collision mo del.
PAGE 80
66 0 5 10 15 20 25 30 35 40 45 20 22 24 26 28 30 32 34 36 R (au)P (au) Figure 4{17. Grid deformation at t = 1400 au, for the dual crossing collision mo del.
PAGE 81
67 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 4 3.5 3 2.5 2 1.5 1 h1 Tlog e (E) (au) EPQCLE SOFFT Tully et al. Figure 4{18. Probabilit y of transmission in the ground state, for the dual crossing collision mo del.
PAGE 82
68 T able 4{3. P arameters used in the NaI complex mo del. P arameter V alue (au) R 0 5.047 P 0 0.0 0.12462 A 1 0.02988 A 2 101. 43 A 12 0.00202 B 2 3.000 C 2 18. 950 + 2.756 12. 179 0.660 E 0 0.07626 1 2.158 12 0.194 R x 13. 24 M 35,482 so that the PWTDM b ecomes, 11 ( P ; R ) = 1 R R 0 2 2 ( P P 0 ) 2 # ; (4.30) with 12 = 21 = 22 = 0. The sim ulation follo ws the resulting motion of this state. The mo del parameters are sho wn in T able 4{3 The diabatic p oten tials are displa y ed in Figure 4{19 4.5.2 Prop erties of In terest In addition to observing the coherence, exp ectation v alue of the p osition and its disp ersion, and the phase space grid, w e also consider b ound and free neutral and ionic p opulations as the NaI complex oscillates from its primarily ionic to primarily co v alen t state. Bound and free p opulations. It is in teresting to follo w the p opulations of the ionic and neutral states as the system ev olv es, giving insigh t in to the nature of the disso ciation in to an asymptotically free neutral system. F or this, w e dene three
PAGE 83
69 0.1 0.05 0 0.05 0 5 10 15 20 25 30 Energy (au)R (au) H 11 H 22 H 12 Figure 4{19. P oten tial curv es for the NaI complex.
PAGE 84
70 p opulations: the ionic, the b ound neutral and the free neutral. The ionic p opulation is the probabilit y of nding the system in the ionic state at an y in tern uclear separation, 2 = 1 Z 1 dR dP 22 ( R ; P ) : (4.31) W e dene the b ound neutral p opulation is the probabilit y of nding the NaI complex in the neutral state at n uclear separation up to the crossing p oin t R = R x b 1 = R x Z 0 dR dP 11 ( R ; P ) ; (4.32) while the free ionic p opulation is dened as the probabilit y of the neutral state from the crossing p oin t b ey ond, f 1 = 1 Z R x dR dP 11 ( R ; P ) : (4.33) This division b et w een b ound and free is motiv ated b y the observ ation that the ma jorit y of the transfer b et w een ionic and neutral p oten tials o ccurs at the a v oided crossing, and that an y part of the w a v epac k et whic h ends up in the neutral state but propagating to w ard innite separation has a negligible probabilit y of shifting to the ionic state m uc h b ey ond the a v oided crossing. 4.5.3 Results The ionic and co v alen t p opulations are displa y ed in Figure 4{20 W e see the oscillations in the p opulations b et w een ionic and co v alen t, rep eating appro ximately ev ery 40 000 au. This pattern can b e compared to the exp ectation v alues of p osition see in Figure 4{21 The p osition oscillates with the same frequency as the c hange in p opulation, sho wing that eac h time the w a v epac k et heads across the a v oided crossing from its b ound co v alen t state, it con v erts almost completely in to the ionic state, with a small amoun t escaping in to the free neutral state. Ov er time, the
PAGE 85
71 free neutral p opulation gradually increases at these crossings, and the NaI slo wly disso ciates. The EPQCLE is quan titativ ely similar to the exact results for the rst half of the sim ulation, and main tains qualitativ e accuracy for the remainder. The disp ersion in Figure 4{21 sho ws an in teresting dierence b et w een the exact and the quan tumclassical algorithm. While the disp ersion con tin ues to rise in the SOFFT sim ulation, it reac hes its rst p eak and then b egins to decline somewhat in the EPQCLE sim ulation. This rerects the nature of the eectiv e p oten tial, where eac h p oin t is guided b y a com bination of excited and ground state forces. Because the ionic curv e do es not p ermit escap e, what w ould normally b e asymptotically free w a v epac k ets tend to b e pulled bac k to w ard the crossing b ecause of the attractiv e ionic p oten tial. Consequen tly the free neutral p opulation is alw a ys lo w er using the EPQCLE equation than the SOFFT, an observ ation supp orted b y Figure 4{20 Since the ma jorit y of the p opulation remains in the ionic or b ound neutral state, and these p opulations are w ell matc hed b et w een the exact and quan tumclassical sim ulations, it is not surprising that the exp ectation v alue of the p osition is quantitativ ely in the b eginning, and qualitativ ely for the remaining, similar for b oth metho ds. Ho w ev er, the disp ersion is m uc h more sensitiv e to the increased asymptotically free w a v epac k ets in the SOFFT sim ulation, and for the reasons discussed, w e nd signican t div ergence b et w een the EPQCLE and SOFFT results. The coherence, sho wn in Figure 4{22 initially p eaks through the rst crossing, but through subsequen t crossings it is substan tially diminished. Ho w ev er, the EPQCLE sho ws quan titativ ely similar results to the SOFFT calculations. The deformation of the phase space grid, plotted in Figure 4{23 has c haracteristics not seen in the other t w o mo dels. One line of p oin ts emits from the cen ter of the cluster, quic kly straigh tening and rerecting the negligible force on the p oin ts. These p oin ts corresp ond to the asymptotically free neutral comp onen ts of the PWTDM. The second group circles around, gaining v elo cit y and p osition, then
PAGE 86
72 turning. These ellipses are c haracteristic of the phase space of classical particles in a w ell, and indeed rerect the quasiclassical motion under the HellmannF eynman force of the PWTDM p oin ts as they follo w the ionic and b ound co v alen t curv es. 4.6 Comparison Using V ariable and Constan t Timesteps In this Section, w e ev aluate the usefulness of the v ariable timestep asp ect of the relaxanddriv e algorithm. T o do this, w e consider the Nasurface algorithm of Section 4.3 and sim ulate using v arying upp er and lo w er tolerances. The n um b er of steps tak en in eac h case is compared with the n um b er that w ould b e required of the same algorithm, but k eeping the timestep xed. The xed timestep w ould necessarily adv ance b y steps no greater than the smallest timestep used in the v ariable timestep algorithm, and it is based on this timestep that w e estimate the corresp onding steps required for the xed timestep approac h. The results are sho wn in Figure 4{24 W e see that as the tolerance decreases (and th us the accuracy increases), the fraction of steps sa v ed b y the in tro duction of the v ariable timestep increases sup erlinearly One concludes that while the v ariable timestep ma y not b e imp ortan t for lo w accuracy sim ulations, considerable computational sa vings can b e had for high accuracy propagation. 4.7 Conclusion By examining three simple t w ostate mo dels, w e w ere able to compare the EPQCLE metho d with the exact SOFFT quan tum mec hanical solution. F or all mo dels, w e found v ery go o d agreemen t b et w een the EPQCLE and the SOFFT results. The agreemen t w as at least qualitativ e, and in man y cases quan titativ e to visual precision. W e also sa w conditions under whic h the quan tumclassical mo del deviated from exact quan tal results. Finally w e compared a xed timestep v arian t of the relaxanddriv e algorithm with the v ariable timestep v ersion. F or the mo del of an alk ali atom approac hing a metal surface, w e examined probabilit y transitions, p osition exp ectation and its deviation, momen tum exp ectation
PAGE 87
73 0 0.5 1 1.5 2 0 20000 40000 60000 80000 100000 120000 Ionic and Neutral Populationstime (au) EPQCLE: Ionic SOFFT: Ionic EPQCLE: Neutral Bound SOFFT: Neutral Bound EPQCLE: Neutral Free SOFFT: Neutral Free Figure 4{20. Ionic and neutral p opulations o v er time, for the NaI complex.
PAGE 88
74 0 5 10 15 20 25 0 10000 20000 30000 40000 50000 60000 70000 80000 90000 100000
PAGE 89
75 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 20000 40000 60000 80000 100000 120000 Coherence: Real( G12)time (au) EPQCLE SOFFT Figure 4{22. Coherence as a function of time, for the NaI complex.
PAGE 90
76 0 20 40 60 80 100 120 140 80 60 40 20 0 20 40 60 80 R (au)P (au) Figure 4{23. Phase space grid at the end of the sim ulation, for the NaI complex.
PAGE 91
77 0 500 1000 1500 2000 2500 3000 3500 10 3 10 1 10 5 10 3 10 7 10 5 10 9 10 7 number of timestepstolerance Variable Fixed Figure 4{24. Num b er of steps required b y the relaxanddriv e algorithm, compared to an estimated n um b er required for a xed timestep v ersion.
PAGE 92
78 and its deviation, densit y function and phase space grid ev olution. W e found that the EPQCLE metho d repro duced quan titativ ely the exact v alues found through the quan tum mo del for all observ ables. W e also found the phase space grid to deform substan tially as the system ev olv ed, initially distorting but ev en tually lining up in a straigh t line as the eectiv e force v anished. F or the dual crossing diabatic surface collision, in addition to probabilit y transitions for a giv en energy w e calculated the transmission probabilit y of the ground state for a wide range of energies. This sho w ed the EPQCLE to deviate at the lo w er energies, where n uclear in terference eects w ere imp ortan t. A t higher energies, the corresp ondence b et w een the mixed quan tumclassical results and the exact quan tum results w as go o d. W e also noted that the phase space grid deformed only for p oin ts whic h sp en t a signican t amoun t of time in the region of strong in teraction. This is reasonable, as the deformation comes from the eectiv e force, whic h is nonv arying (and in fact zero) outside the in teraction. The NaI mo del sho w ed in teresting bac kandforth transitions b et w een the ionic and co v alen t state, as the in tern uclear distance oscillated due to the longrange attractiv e ionic p oten tial. W e also sa w that eac h passing through the a v oided crossing led to a small amoun t of disso ciation in to the asymptotically free neutral state. Ionic, b ound neutral and free neutral p opulations w ere quan titativ ely similar for the EPQCLE and SOFFT algorithms for the b eginning of the sim ulations, and qualitativ ely similar for the remainder. Coherences w ere quan titativ ely similar, sho wing an initial p eak after the rst pass through the a v oided crossing, and m uc h smaller magnitudes thereafter. Finally the grid deformation sho w ed the grid to split in to t w o groups. The rst represen ted the asymptotically free neutral state, where the p oin ts form a straigh t line as the HellmannF eynman force v anishes and the particles propagate to innite distances. The second group circled in phase space, as the
PAGE 93
79 p oin ts oscillated bac k and forth in the p oten tial w ell formed b y the ionic attraction at large distances and the neutral repulsion at small distances. The EPQCLE metho d w as sho wn to b e v ery robust under a wide v ariet y of conditions. Comparing its v ariable timestep v ersion (whic h w as used in all the sim ulations) to the xed timestep alternativ e, w e found that not only w as the EPQCLE metho d accurate, but use of the relaxanddriv e propagation generated m uc h more ecien t sim ulations than could ha v e b een obtained without dynamically v arying the timestep.
PAGE 94
CHAPTER 5 ALKALI A TOMRARE GAS CLUSTERS: GENERAL F ORMULA TION 5.1 In tro duction In this c hapter, w e describ e the in teractions of an alk ali atom (Ak) em b edded in a cluster of rare gas atoms (Rg N ). Cluster dynamics are in teresting insofar as they presen t a bridge b et w een isolated atoms and the (eectiv ely) innite bulk liquid, and can pro vide insigh t in to the dynamics of molecules on extended surfaces. The AkRg N clusters are among the simplest cluster systems to study b ecause of the presence of a single v alence electron in the alk ali atom, and the closed shell structure of the rare gas atoms. Pseudop oten tial descriptions of the in teractions b et w een the v alence electron ( e ), the Ak core and the Rg atoms can greatly simplify the description of the cluster with v ery mo dest p enalties in accuracy In the con text of mixed quan tumclassical mo dels, these clusters are w ell suited to sim ulation through the EPQCLE b y treating the v alence electron quan tally and the n uclear cores quasiclassically F or the ma jorit y of the c hapter, w e consider the general case of an alk ali atom em b edded in a rare gas cluster, although our sim ulations in Chapter 6 fo cus on the sp ecic case of the lithium atom (Li) em b edded in a helium cluster (He N ). 5.2 Ph ysical System W e consider a cluster initially at thermal equilibrium, and in tro duce a ground or excited alk ali atom to the cen ter of the cluster. W e do not concern ourselv es with the means b y whic h the alk ali atom is excited or em b edded within the cluster, although t ypically the excitation is due to a laser pulse. F rom this initial setup, w e follo w the dynamics of the alk ali atom as it deexcites and mo v es within the cluster. W e further assume the cluster is isolated from an y en vironmen tal eects, so that sp ectral and conguration measuremen ts represen t 80
PAGE 95
81 those of an isolated AkRg N cluster. As w e see later, w e in tro duce a con taining p oten tial to k eep the Rg atoms from disp ersing. This p oten tial is used strictly to main tain a cluster formation, and do es not represen t in teraction with an en vironmen t. 5.3 Prop erties of In terest First of all, w e are in terested in the structure of the Rg cluster in equilibrium, b efore the in tro duction of the Ak atom. W e need to repro duce kno wn densit y proles and pairpair correlation functions to ensure the cluster is represen tativ e of a real ph ysical system. Secondly w e are in terested in follo wing the dynamics of the Ak atom as it mo v es within the cluster. In particular, w e are lo oking for migration from the cen ter of the cluster to its surface, as other exp erimen tal and theoretical studies indicate that the Ak atom tends to reside on the cluster surface. Finally w e wish to compute the time dep enden t emission sp ectra resulting from the deexcitation of the Ak atom, Ak ( n 0 l 0 ) + R g N Ak ( nl ) + R g N + ; (5.1) where is a photon with energy corresp onding to the electronic deca y 5.4 Hamiltonian for Alk aliRare Gas P airs F ollo wing the dynamics of an AkRg N cluster through a full quan tal treatmen t w ould b e an extremely demanding problem computationally for more than a few Rg atoms. Instead, w e reduce the AkRg in teraction to a threeb o dy problem b y using pseudop oten tial in teractions b et w een the n uclear cores and the electron. The pseudop oten tial treatmen t is explored extensiv ely for the single AkRg pair in Rey es. 109 In this section, w e summarize this approac h. In order to describ e the threeb o dy in teraction, w e consider a xed n uclear conguration with R AB the p osition v ector from the alk ali core (A) to the rare gas atom (B), and r A ( r B ) the p osition v ector from the alk ali core (rare gas atom) to
PAGE 96
82 the electron e With this notation, w e can write the Hamiltonian in v e distinct comp onen ts, ^ H pair = 1 2 r 2r A + ^ V Ak A ( r A ) + ^ V Rg B ( r B ) + ^ V cr oss AB ( r A ; R AB ) + ^ V cor e AB ( R AB ) : (5.2) W e ha v e implicitly tak en the PWT o v er n uclear v ariables but not electronic v ariables. Because the p oten tials can all b e expressed as p olynomial functions of ^ Q the PWT amoun ts to replacing ^ Q b y R throughout the Hamiltonian. In what follo ws, w e will assume w e are w orking with partially Wigner transformed op erators, but drop the subscript 'W' for notational simplicit y Also note that the electronic v ariable remains quan tal, although w e will use the notation r rather than ^ q for consistency with common usage in the literature. The rst term on the RHS of Eq. 5.2 is, of course, the kinetic energy op erator of the v alence electron. The second and third terms are the p oten tials arising from the in teraction of the electron with the Ak core and the Rg atom, resp ectiv ely The ^ V cr oss AB term is a crossterm stemming from the p olarization of the Rg atom b y b oth the v alence electron and Ak core. Finally the last term ^ V cor e AB is the in teraction b et w een the Ak core and the Rg atom. W e examine eac h of these p oten tials in detail for the general AkRg pair, and pro vide sp ecic parameters for the LiHe in teraction in Chapter 6. The e Ak core p oten tial can b e divided in to three comp onen ts, ^ V Ak A ( r A ) = Z A r A + ^ V pol Ak A ( r A ) + ^ V sr A ( r A ) ; (5.3) where the rst term is the Coulom b in teraction b et w een the Ak core of c harge Z A and the electron. The second con tribution arises from the dip ole p olarization of the
PAGE 97
83 Ak core b y the v alence electron, ^ V pol Ak A ( r A ) = d A 2 r 4 A w ( r A ; A ) 2 ; (5.4) with d A the dip ole p olarizabilit y of the core and w a cuto function of the distance w ( r ; ) = 1 exp ( r 2 ) : (5.5) The nal con tribution to ^ V Ak A is an l dep enden t shortrange pseudop oten tial, ^ V sr X ( r X ) = X l ;i B l ;i exp ( l ;i r 2 X ) P l ;X ; (5.6) where B l ;i and l ;i are pseudop oten tial parameters adjusted to t exp erimen tal data, P l ;X is the pro jection op erator on angular symmetry l P l ;X = X m j Y l m ( ^ r X ) ih Y l m ( ^ r X ) j ; (5.7) and the states fj Y l m ( ^ r X ) ig are the spherical harmonic functions cen tered on the core X. This pseudop oten tial sim ulates repulsion due to the eects of the P auli exclusion principle when the v alence electron approac hes the core electrons, as w ell as the attraction due to the incomplete screening of the n uclear c harge. The third term in Eq. 5.2 is the e Rg p oten tial, and can b e divided in to t w o terms, ^ V Rg B ( r B ) = ^ V pol Rg B ( r B ) + ^ V sr B ( r B ) : (5.8) The rst comp onen t stems from the p olarization of the Rg atom b y the electron, ^ V pol Rg B ( r B ) = d B 2 r 4 B w ( r B ; ) 4 q 0 B 2 r 6 B w ( r B ; ) 6 ; (5.9)
PAGE 98
84 where q 0 B = q B 6 1 q B b eing the quadrup ole p olarizabilit y of the Rg atom and 6 1 the dynamical correction to the static p olarizabilit y The shortrange pseudop oten tial ^ V sr B is dened in Eq. 5.6 The fourth term in Eq. 5.2 is a crossterm arising from the p olarization of the Rg atom b y b oth the Ak core and the v alence electron, ^ V cr oss AB ( r A ; R AB ) = Z A d B P 1 ( cos ) R 2 AB r 2 B w ( r B ; ) 2 + Z A q B P 2 ( cos ) R 2 AB r 2 B w ( r B ; ) 3 : (5.10) Here, P 1 and P 2 are the Legendre p olynomials and is the angle b et w een the v ectors R AB and r B Finally the last term in Eq. 5.2 is the in teraction b et w een the alk ali core and the rare gas atom. It is assumed to in v olv e shortrange, dip ole and quadrup ole con tributions, ^ V cor e AB ( R AB ) = ^ V sr AB ( R AB ) 1 2 d B ( R 2 AB + d 2B ) 2 1 2 q 00 B ( R 2 AB + d 2B ) 3 ; (5.11) where q 00 B = q B + 2 d B d 2B The shortrange term ^ V sr AB has b een determined b y assuming it to ha v e the form, ^ V sr AB ( R AB ) = A exp ( bR AB ) ; (5.12) and adjusting the parameters A and b so that ^ V cor e AB ( R AB ) in Eq. 5.11 ts the most accurate curv es in the literature. This pro cedure has b een sho wn to accurately repro duce AkRg energy curv es for a v ariet y of alk ali atom and rare gas com binations. 109 5.5 Hamiltonian for the Alk aliRare Gas Cluster The cluster is similar to the pair, in that eac h AkRg pair in the cluster is treated as a threeb o dy problem in v olving alk ali core, v alence electron and rare gas atom, and describ ed using the same Hamiltonian elemen ts. In the case of the cluster, ho w ev er, w e also require that the Rg atoms main tain a cohesiv e structure
PAGE 99
85 rather than drift apart. T o ensure this, w e imp ose a constraining p oten tial on the system. Dening R AB i as the p osition v ector from the alk ali core to the Rg atom i and R B i the p osition v ector of the Rg atom from the origin, w e can write the cluster Hamiltonian as, ^ H cl uster = 1 2 r 2r A + ^ V Ak A ( r A ) + N X i =1 ^ V Rg B i ( r B ) + N X i =1 h ^ V cr oss AB i ( r A ; R AB i ) + ^ V cor e AB i ( R AB i ) + ^ V hol d i ( R B i ) i + N X j >i =1 ^ V cor e ( j R B i R B j j ) : (5.13) The nal term in Eq. 5.13 is the in teraction b et w een the Rg atoms for the quasiclassical motion. F ollo wing Aziz and co w ork ers, 110 w e tak e the RgRg p oten tial to ha v e the form, for in tern uclear distance R ^ V cor e ( R ) = V ( x ) ; (5.14) where V ( x ) = A exp ( x + x 2 ) F ( x ) X j =0 2 C 2 j +6 x 2 j +6 ; (5.15) F ( x ) = 8><>: exp [ ( D =x 1) 2 ] ; x D 1 ; x > D : (5.16) This p oten tial has b een written in terms of the dimensionless distance x = R =R M The constraining p oten tial, ^ V hol d i is tak en to b e a sigmoidal function cen tered along the b oundary of a sphere of radius R hol d ^ V hol d i ( R B i ) = a 1 + exp [ b ( R B i R hol d )] : (5.17)
PAGE 100
86 This p oten tial sho ws asymptotes ^ V hol d i ( 1 ) = 0 and ^ V hol d i (+ 1 ) = a with midp oin t at the holding radius ^ V hol d i ( R hol d ) = a= 2. The steepness of the function is con trolled b y the parameter b and determines the strength and range of the holding force. It acts on the Rg atoms only k eeping them b ound roughly within a sphere of radius R hol d while p ermitting the alk ali atom free motion within the cluster. 5.6 Electronic Sp ectral Calculations When a molecular system undergo es electronic motion, the accelerating c harges emit electromagnetic radiation. A t distances large compared to the electronic motion, the rux of this radiation at p oin t r from the cen ter of the system is giv en b y the P o yn ting v ector, 111 S ( t ) = 1 0 ( E B ) = 0 16 2 c [ D ( t )] 2 sin 2 r 2 ^ r ; (5.18) where D ( t ) is the dip ole momen t of the system, 0 is the p ermittivit y constan t and c is the sp eed of ligh t. By in tegrating o v er all angles, w e obtain the p o w er emited b y the source, P ( t ) = Z n S ( t ) d a = j A ( t ) j 2 (5.19) where A ( t ) = r 1 6 0 c 3 D ( t ) : (5.20) F rom Eq. 5.20 w e can see that the computation of the dip ole momen t is essen tial to the sp ectral calculation. In a mixed quan tumclassical system, the dip ole momen t is obtained b y calculating the exp ectation v alue of the dip ole op erator ^ D D ( t ) = T r [ ^ ( t ) ^ D ] = T r[( t )] : (5.21)
PAGE 101
87 F or the AkRg N cluster, the dip ole op erator can b e written with the origin at the Ak core, so that it b ecomes, ^ D = r + d A r A r 3A w ( r A ; A ) + d B N X i =1 r B i r 3B i w ( r B i ; B ) + Z A R AB R 3AB w ( R AB ; B ) # : (5.22) The rst term is the dip ole from the v alence electron. The second is the induced dip ole in the Ak core b y the v alence electron. The third and fourth terms sum o v er all Rg atoms, and represen t the induced dip ole in the Rg core b y the v alence electron and Ak core, resp ectiv ely In practice, the induced dip oles are m uc h smaller than the v alence electron dip ole. W e can calculate the emission sp ectrum b y taking the F ourier transform of Eq. 5.19 I ( ) = 1 p 2 Z P ( t ) exp( i! t ) dt = 1 p 2 Z j A ( t ) j 2 exp ( i! t ) dt: (5.23) If the emission sp ectrum c hanges o v er time, one can use a windo w ed F ourier transform cen tered ab out the time of in terest. In the case of AkRg N clusters, this allo ws us to trace the ev olution of the emission sp ectrum as the alk ali atom mo v es from the cen ter of the cluster to its surface. 5.7 Electronic Basis of Gaussian A tomic F unctions 5.7.1 Equations of Motion W e recall the EPQCLE in an arbitrary basis, Eq. 3.7 repro duced here for con v enience: d dt = ( i H q n y ) S 1 S 1 ( i H q + n ) ; (5.24)
PAGE 102
88 where w e ha v e used the notation, n h j @ =@ t j i + dR dt h j @ =@ R j i + dP dt h j @ =@ P j i : (5.25) F or the cluster, H q = H cl uster A con v enien t basis in v olv es Gaussian functions cen tered on the Ak core, and leads to matrix elemen ts whic h can b e ev aluated analytically T o this end, w e in tro duce primitiv e cartesian Gaussian states, dened b y an in teger triplet a = ( a x ; a y ; a z ), a Gaussian cen ter A = ( A x ; A y ; A z ) and an exp onen tial co ecien t a suc h that their pro jection on the electronic co ordinate r is h r j a ; A ; a i = ( r x A x ) a x ( r y A y ) a y ( r z A z ) a z exp [ a ( r A ) 2 ] : (5.26) F rom these primitiv es, w e can form Gaussian atomic functions (GAFs) whic h resemble the h ydrogenic orbitals. The lo w est three symmetries are j s ; R ; i j a = (0 ; 0 ; 0); R ; i = j 000; R ; i ; (5.27) j p ; R ; i 2 fj 100; R ; i ; j 010; R i ; j 001; R ; ig ; (5.28) j d ; R ; i 2 fj 110; R ; i ; j 011; R ; i ; j 101 ; R ; i ; j 200 ; R ; i j 020; R ; i ; 2 j 002; R ; i j 200; R ; i j 020; R ; ig : (5.29) Finally w e construct the basis elemen ts b y com bining linear com binations of these GAFs in segmen ted con tractions, j i = N X j c j j ; R ; j i ; (5.30) where 2 f s; p; d g f c j g are the con traction co ecien ts and N is a normalization factor. W e can write the full M dimensional basis as a ro w v ector of these segmen ted
PAGE 103
89 con tractions, j i = ( j 1 i j 2 i : : : j M i ) ; (5.31) corresp onding to the notation used to deriv e Eq. 3.7 Since j i do es not dep end explicitly on time or n uclear momen ta, n reduces to a momen tum coupling term, n = dR dt h j @ @ R j i : (5.32) A non trivial matter is the analytical computation of these matrix elemen ts, atten tion to whic h w e turn presen tly 5.7.2 Ov erlap Matrix Elemen ts In computing matrix elemen ts, w e fo cus on the essen tial part of the computation, whic h is the matrix elemen t b et w een primitiv e cartesian Gaussians. F or the o v erlap S w e need to compute h a j b i whic h is p ossible using the recursion form ulae from Obara, 112 113 h a + 1 i j b i = ( P i A i ) h a j b i + 1 2 a i h a 1 i j b i + 1 2 b i h a j b 1 i i ; (5.33) h a j b + 1 i i = ( P i B i ) h a j b i + 1 2 h a j b 1 i i + 1 2 a i h a 1 i j b i ; (5.34) whhere 1 i is an in teger triplet ( ix ; iy ; iz ), = a + b and P = ( a A + b B ) = Recursion con tin ues to the base case, h s ; A ; a j s ; B ; b i = 2 exp [ ~ ( A B ) 2 ] ; (5.35) where ~ = a b = If computational eciency is a concern, some of the lo w er recursion relations can b e rolled out explicitly rather than computed recursiv ely
PAGE 104
90 5.7.3 Kinetic Energy Matrix Elemen ts F or the kinetic energy w e wish to compute h a j ^ K j b i = h a j 1 2 r 2r j b i W e can do this b y using the recursion relations, h a + 1 i j ^ K j b i = ( P i A i ) h a j ^ K j b i + 1 2 a i h a 1 i j ^ K j b i + 1 2 b i h a j ^ K j b 1 i i = 2 ~ h a + 1 i j b i 1 2 a a i h a 1 i j b i ; (5.36) h a j ^ K j b + 1 i i = ( P i B i ) h a j ^ K j b i + 1 2 b i h a j ^ K j b 1 i i + 1 2 a i h a 1 i j ^ K j b i = 2 ~ h a j b + 1 i i 1 2 b b i h a j b 1 i i : (5.37) The recursion con tin ues to the base case, h s ; A ; a j ^ K j s ; B ; b i = ~ (3 2 ~ ( A B ) 2 ) h s ; A ; a j s ; B ; b i = ~ (3 2 ~ ( A B ) 2 ) 2 exp [ ~ ( A B ) 2 ] : (5.38) 5.7.4 Coulom b Matrix Elemen ts The n uclear attraction in tegral is a threecen ter in tegral, h a j ^ V C j b i = Z d r ( A ; a ) 1 j r C j ( B ; b ) : (5.39) F ollo wing Obara, 113 this can b e ev aluated using the recursion relations, h a + 1 i j ^ V C j b i ( m ) = ( P i A i ) h a j ^ V C j b i ( m ) ( P i C i ) h a j ^ V C j b i ( m +1) + 1 2 a i h h a 1 i j ^ V C j b i ( m ) h a 1 i j ^ V C j b i ( m +1) i + 1 2 b i h h a j ^ V C b 1 i i ( m ) hj ^ V C j b 1 i i ( m +1) i ; (5.40) and similarly for h a j ^ V C j b + 1 i i The recursion base is found to b e, h s ; A ; a j ^ V C j s ; B ; b i ( m ) = 2 ~ 1 = 2 h s ; A ; a j s ; B ; b i F m ( U ) ; (5.41)
PAGE 105
91 where U = ( P C ) 2 and F m is the sp ecial function, F m ( T ) = 1 Z 0 dt t 2 m exp ( T t 2 ) : (5.42) 5.7.5 Momen tum Coupling Matrix Elemen ts The momen tum coupling elemen ts h j @ =@ R j i are ev aluated in terms of the o v erlap elemen ts, h a j @ =@ B x j b i = b h a j b + 1 x i b 1 h a j b 1 x i ; (5.43) h a j @ =@ B y j b i = b h a j b + 1 y i b 2 h a j b 1 y i ; (5.44) h a j @ =@ B z j b i = b h a j b + 1 z i b 3 h a j b 1 z i : (5.45) Eac h elemen t on the RHS of Eqs. 5.43 to 5.45 is a straigh tforw ard o v erlap in tegral, whose ev aluation w e ha v e already deriv ed previously 5.7.6 Dip ole Matrix Elemen ts Recalling the dip ole for the AkRg system, w e see that w e require the follo wing matrix elemen ts, h a j r A j b i ; (5.46) h a j r A r 3A (1 exp [ r A A ]) j b i ; (5.47) h a j r B r 3B (1 exp [ r B B ]) j b i ; (5.48) h a j b i : (5.49) The rst in tegral is a momen t in tegral, whic h can b e computed using the recursiv e pro cedure from Obara. 113 Consider the general momen t in tegral, h a jM ( ) j b i ; (5.50)
PAGE 106
92 where M ( ) = x x y y z z : (5.51) These in tegrals can b e found through the recursiv e form ula, h a + 1 i jM ( ) j b ) = ( P i A i ) h a jM ( ) j b i + 1 2 a i h a 1 i jM ( ) j b i + 1 2 b i h a jM ( ) j b 1 i i + 1 2 i h a jM ( 1 i ) j b i : (5.52) The recursion con tin ues to the base case, h s ; A ; a jM ( ~ 0 ) j s ; B ; b i = h s ; A ; a j s ; B ; b i = 2 exp [ ~ ( A B ) 2 ] : (5.53) The second and third in tegrals can b e calculated using a pro cedure presen ted b y McMurc hie and Da vidson, 114 and summarized in detail in Rey es. 109 The nal in tegral is a simple o v erlap in tegral, b ecause the op erator do es not dep end on electronic co ordinates. W e ha v e presen ted the calculation of the o v erlap in tegral previously 5.7.7 Pseudop oten tial Matrix Elemen ts In order to calculate the pseudop oten tial matrix, w e need the additional elemen ts, h a j 1 r 4 [1 exp ( c r 2 )] 2 j b i ; (5.54) h a j 1 r 4 [1 exp ( c r 2 )] 4 j b i ; (5.55) h a j 1 r 6 [1 exp ( c r 2 )] 6 j b i ; (5.56) h a j l X m = l j l m ih l m j j b i : (5.57)
PAGE 107
93 The rst three elemen ts are computed using a pro cedure from McMurc hie and Da vidson, 114 while the fourth elemen t is calculated based on metho ds in Sc h w erdtfeger and Silb erbac h. 115 These calculations are summarized in detail in Rey es. 109 5.8 Computing the Quasiclassical T ra jectory The c hange in p osition of core i (Ak or Rg) is simply calculated from the momen tum, d R i dt = P i M i : (5.58) Ho w ev er, the c hange in momen tum is substan tially more in v olv ed. Eac h quasiclassical v ariable k is sub jected to the HellmannF eynman force, F k = T r ^ @ ^ H @ R k # = T r[] ; (5.59) where w e ha v e used the abbreviated notation ^ H = ^ H cl uster Using our nonorthogonal basis j i this expression b ecomes, F k = T r h j @ ^ H @ R k j i # = T r[ S ] = 1 T r [ S ] T r[ @ H @ R k ] T r[ n yR k S 1 H ] T r [ HS 1 n R k ] : (5.60) The rst term on the RHS of Eq. 5.60 in v olv es the partial deriv ativ e of the cluster Hamiltonian matrix, and can b e ev aluated n umerically b y computing H at a p erturb ed co ordinate R k + @ H @ R k = H ( R k + ) H ( R k ) : (5.61) 5.9 Computational Details Examining the EPQCLE, w e see that a n um b er of matrices can b e computed a single time at the b eginning of the sim ulation, suc h as the o v erlap, momen tum, electronic dip ole and kinetic energy Ho w ev er, the pseudop oten tial is dep enden t
PAGE 108
94 on the relativ e AkRg congurations, and m ust b e computed at eac h timestep. In addition, the n umerical ev aluation of the HellmannF eynman force requires that the Hamiltonian for sev eral p erturb ed congurations b e ev aluated for eac h force calculation. Consequen tly it is imp ortan t that the pseudop oten tial calculations run quic kly This can b e done b y precomputing the pseudop oten tial matrix for sev eral in teratomic distances along the zaxis, extrap olating to the desired distance, and then rotating to the actual orien tation. In this section, w e demonstrate that the rotation is a simple matter of m ultiplication b y rotation matrices. Consider a singlecen ter basis j i and general op erator ^ O in the threeb o dy AkRg system, where the Ak atom is cen tered at the origin, R is the p osition v ector from this core to the Rg atom, and r is the p osition v ector of the v alence electron to the Ak core. Supp ose that w e ha v e precomputed O = h j ^ O j i for sev eral v alues of R along the zaxis. Our task is to compute the new matrix O 0 for the general case where the Rg atom has b een rotated b y n = ( ; ), where and are the p olar and azim uthal angles, resp ectiv ely If the op erator is of the form, ^ O = ^ O ( j R j ; j r j ; j R r j ) ; (5.62) whic h is the case for the cluster Hamiltonian op erator, w e can compute O 0 through a rotation D (n), j 0 i = j i D (n) : (5.63) Algebraically w e can see this b y examining a single matrix elemen t with R along the zaxis, O = h j ^ O j i = Z d r ( r ) ^ O ( j R j ; j r j ; j R r j ) ( r ) : (5.64)
PAGE 109
95 When R is in the rotated p osition R 0 = D (n) R O 0 = Z dr ( r ) ^ O ( jD (n) R j ; j r j ; jD (n) R r j ) ( r ) : (5.65) Substituting r D (n) r = D r O 0 = Z dr ( D r ) ^ O ( jD R j ; jD r j ; jD R D r j ) ( D r ) : (5.66) Since a rotation of a v ector do es not c hange its magnitude, O 0 = Z dr ( D r ) ^ O ( j R j ; j r j ; j R r j ) ( D r ) (5.67) ) O 0 = Z dr ( D r ) ^ O ( j R j ; j r j ; j R r j ) ( D r ) = D y h j ^ O j i D : (5.68) T urning our atten tion to the construction of D w e use the fact that our basis in v olv es segmen ted con tractions of Gaussian atomic functions. In this case, basis elemen ts of same symmetry mix, while dieren t symmetries do not. W e will sho w the full deriv ation for s t yp es and p t yp es. Denoting the s t yp e con traction, s ( r ) = X i exp ( i r 2 ) ; (5.69) w e w an t to calculate s ( D r ). Since D r implies the follo wing co ordinate transformations, x x cos cos + y sin sin z sin ; (5.70) y x sin + y cos ; (5.71) z x cos sin + y sin sin + z cos : (5.72)
PAGE 110
96 Then, s ( D r ) = X i c i exp ( i r 2 ) = s ( r ) : (5.73) Similarly consider the p t yp e con tractions, p x ( r ) = X i c i x exp ( i r 2 ) ; (5.74) p y ( r ) = X i c i y exp ( i r 2 ) ; (5.75) p z ( r ) = X i c i z exp ( i r 2 ) : (5.76) Then p x ( D r ) = X i c i ( x cos cos + y sin cos z sin ) exp( i r 2 ) = p x ( r ) cos cos + p y ( r ) sin cos + p z ( r )( sin ) ; (5.77) p y ( D r ) = p x ( r )( sin ) + p y ( r ) cos ; (5.78) p z ( D r ) = p x ( r ) cos sin + p y ( r ) sin sin + p z ( r ) cos : (5.79) Supp ose our basis is ( r ) = [ s ( r ) p x ( r ) p y ( r ) p z ( r )] : (5.80) Then, ( D r ) = [ s ( D r ) p x ( D r ) p y ( D r ) p z ( D r )] = [ s ( r ) p x ( r ) p y ( r ) p z ( r )] D (n) = ( r ) D sp (n) ; (5.81)
PAGE 111
97 where D sp (n) = 0BBBBBBB@ 1 0 0 0 0 cos cos sin cos sin 0 sin cos cos sin sin 0 sin 0 cos 1CCCCCCCA : (5.82) W e can w ork out the mixing of d functions similarly the results of whic h w e presen t in T able 5.9 5.10 Conclusion W e ha v e thoroughly describ ed the metho ds to ev olv e the AkRg N cluster, starting with the Hamiltonian for the AkRg pair and generalizing to a full cluster in three dimensions. W e ha v e treated the Ak v alence electron explicitly and ha v e computed the electron in teractions using semilo cal l dep enden t pseudop oten tials. In the context of the EPQCLE w e ha v e treated the electron as a quan tal v ariable, and the Ak and Rg cores as quasiclassical. By computing the dip ole of the molecular system, w e ha v e sho wn ho w the electronic emission sp ectra of an initially excited Ak atom can b e computed. In order to carry out the ev olution of the system n umerically w e ha v e in tro duced a basis of segmen ted con tractions of Gaussian atomic functions, and ha v e ev aluated all required matrix elemen ts explicitly The n umerical computation of the HellmannF eynman force guiding the quasiclassical motion of the n uclear cores has b een deriv ed, and the greatly accelerated computation of the pseudop oten tial matrices through table lo okup, n umerical in terp olation and rotation has b een discussed in detail.
PAGE 112
98 T able 5{1. Pseudop oten tial rotation for d function mixing. D ij (n) elemen t v alue D xy ;xy (2 cos 2 1) cos D xy ;y z (2 cos 2 1) sin D xy ;xz 2 sin cos sin cos D xy ;x 2 y 2 2 sin cos cos 2 + sin cos sin 2 D xy ; 2 z 2 x 2 y 2 p 3 sin cos sin 2 D y z ;xy cos sin D y z ;y z cos cos D y z ;xz 2(cos 2 1) sin D y z ;x 2 y 2 sin sin cos D y z ; 2 z 2 x 2 y 2 p 3 sin sin cos D xz ;xy sin sin D xz ;y z sin cos D xz ;xz (2 cos 2 1) cos D xz ;x 2 y 2 cos sin cos D xz ; 2 z 2 x 2 y 2 p 3 cos sin cos D x 2 y 2 ;xy 2 sin cos cos D x 2 y 2 ;y z 2 sin cos sin D x 2 y 2 ;xz (2 cos 2 1) sin cos D x 2 y 2 ;x 2 y 2 (2 cos 2 1) cos 2 + 1 2 (2 cos 2 1) sin 2 D x 2 y 2 ; 2 z 2 x 2 y 2 q 3 4 (2 cos 2 1) sin 2 D 2 z 2 x 2 y 2 ;xy 0 D 2 z 2 x 2 y 2 ;y z 0 D 2 z 2 x 2 y 2 ;xz p 3 sin cos D 2 z 2 x 2 y 2 ;x 2 y 2 q 3 4 sin 2 D 2 z 2 x 2 y 2 ; 2 z 2 x 2 y 2 cos 2 1 2 sin 2
PAGE 113
CHAPTER 6 LITHIUMHELIUM CLUSTERS 6.1 In tro duction In Chapter 5, w e explored all asp ects of the sim ulation of alk alirare gas clusters, in the con text of the EPQCLE. In this c hapter, w e rene our description to the sp ecic case of a lithium atom em b edded in a cluster of helium atoms. The formalism is the same, in that the helium atoms and lithium core are treated quasiclassically while the lithium v alence electron is describ ed quan tally These clusters are only stable at v ery lo w temp eratures (a few degrees Kelvin), where core quan tal eects b ecome imp ortan t. Ho w ev er, b y appro ximating the heliumhelium in teraction with an eectiv e p oten tial, w e are able to mo del bulk liquid helium classically suc h that its temp erature, densit y and radial distribution function matc h exp erimen tal and path in tegral results. F urthermore, b y extracting a sphere of helium atoms from this bulk, and imp osing an appropriate constraining p oten tial, w e are able to repro duce a liquid helium droplet whose densit y prole matc hes quan tum Mon te Carlo sim ulations. With an adequate classical description of the cluster, w e in tro duce a ground state lithium atom in to its cen ter and monitor the conguration and electronic p opulation dynamics o v er time. W e also study the ev olution of the excited lithium atom on a mo del of the cluster surface, in tro ducing electromagnetic elds to induce electronic transitions and ligh t emission. 6.2 Description of the System There are t w o stable isotop es of liquid helium, fermionic 3 He and b osonic 4 He. 116 W e are fo cused on the 4 He isotop e (whic h w e will designate as simply He), although in the classical appro ximation the only dierence resides in their mass. A t temp eratures b elo w 4.2 K, helium v ap or condenses in to the liquid state, 117 while b elo w 99
PAGE 114
100 2.17 K the liquid exhibits sup erruid prop erties. 118 117 Indeed, the thermal w a v elength at these energies exceeds the t ypical corecore separation, and quan tal eects suc h as zerop oin t motion b ecome imp ortan t. 119 Ho w ev er, the excited electronic state of the He atom has an energy of 2 : 3 10 5 K ab o v e the ground state, making it v ery reasonable to assume the He atoms remain in their ground electronic conguration. 118 A t these temp eratures, the lo w kinetic energy p ermits the formation of stable clusters of helium con taining sev eral tens to sev eral h undreds of He atoms. These liquid droplets exhibit the bulk densit y and structure near their cen ter, with a decreasing densit y to w ard their surface. 120 When lithium atoms are in tro duced in to the cluster, they tend to reside on the surface, b ound b y a v ery shallo w w ell in the ground electronic state. 13 When the surface Li atom is excited, its subsequen t b eha vior dep ends hea vily on the excited state. The dierence in excited states is pictorally describ ed in Figure 6{1 Near a rat helium surface, the rst t w o excited states are degenerate, and rerect the p orbital aligned parallel to the surface. This Li(2 p ) conguration minimizes the electronic o v erlap of the helium with the lithium, and results in an attraction to w ard the surface. The third excited state, 2 p has the p orbital aligned orthogonally to the surface, where the electronic o v erlap induces repulsion b et w een the atoms. 121 While the Li(2 p ) tends to mo v e a w a y from the surface with minimal distortion of the He distribution, the Li(2 p ) mo v es to w ard the surface, where the He atoms resp ond b y clustering around the attractiv e Li atom. In exp erimen tal studies with excited Na, excimers of Na and He are found to desorb from the He surface within 70 to 700 ps, dep ending on the initial attractiv e state. 121 The desorb ed excimer is then found to emit to the red of the gas phase Na(3 p 3 s ) transition. Indeed, Li(2 p ) surrounded b y surface He atoms also emits to the red of the corresp onding
PAGE 115
101 A B Figure 6{1. Sc hematic of Li(2 p ) ab o v e a He surface. A) Li(2 p ). B) Li(2 p ). gas phase transition, but in terestingly no ligh t emission is observ ed when the Li(2 s ) is excited in the He bulk. 122 6.3 Prop erties to b e In v estigated First of all, w e are in terested in a classical description of the He atoms. While quan tal n uclear eects are clearly imp ortan t for ab initio descriptions of liquid helium, w e are fo cused on the in teraction of Li with the He, and seek only to reproduce the quan tum exp ectation of the He conguration. T o this end, it is imp ortan t to nd an eectiv e p oten tial for the HeHe in teraction that generates an appro ximately correct He distribution through classical sim ulation, for a giv en He densit y and temp erature. Second, w e wish to mo del the liquid He droplet, where the He atoms hold together and exhibit a densit y distribution that approac hes the bulk helium densit y in the cen ter and tap ers o to w ard the edge. This can b e done b y using an appropriate constraining p oten tial whic h k eeps the atoms in a cluster. Ha ving formed a classical mo del of the He droplet, w e wish to examine the dynamics of a Li atom in tro duced in to the cluster. Sp ecically w e are in terested in
PAGE 116
102 follo wing the n uclear motion and electronic p opulations of ground state Li initially em b edded in the cen ter of the cluster, to b etter understand the mec hanism b y whic h Li atoms tend to reside on the surface of the He cluster. Once w e ha v e clearly established that the Li atom is not miscible with the He liquid, w e shift fo cus to the in teraction b et w een Li and the surface He atoms. In part, w e are in terested in the congurational ev olution of the He atoms and Li core once the Li atom is elev ated to an excited electronic state. Dep ending on whether the or state is p opulated, the dynamics are exp ected to b e quite dieren t, and these dierences are in v estigated. W e are also in terested in the c hange in electronic p opulations and energy lev els as the excited Li atom either approac hes or recedes from the He surface. By in troducing a classical electromagnetic eld resonan t to electronic transitions, w e are able to induce p opulation transitions in Li(2 p ), and stim ulate ligh t emission that can b e seen in the dip ole emission sp ectrum. In addition to pro viding sp ectral results for the LiHe N cluster, w e gain in teresting insigh ts in to the inclusion of an external electromagnetic eld within the mixed quan tumclassical con text. 6.4 Preparation of LithiumHelium Clusters 6.4.1 Bulk Helium Our rst task is to generate a reasonable classical sim ulation of liquid He at ultralo w temp eratures. The four principal elemen ts in v olv ed in this kind of sim ulation are the system b oundaries, the propagation sc heme, the equilibration to attain the correct temp erature, and the in teratomic p oten tials. W e will discuss eac h of these in turn. In order to capture the eectiv ely innite spatial exten t of the liquid on the atomic scale, an attractiv e approac h is to use p erio dic b oundary conditions. 34 In this sc heme, a cen tral b o x is (imagined to b e) replicated at its sides and corners, so that there are 26 additional b o xes surrounding the cen tral v olume. Eac h of these
PAGE 117
103 b o xes con tains precisely the same conguration of He atoms as the cen tral b o x, and an y mo v emen t of an atom in the cen tral b o x (considered a r e al atom) is replicated b y the atoms in the surrounding b o xes (considered virtual atoms). If the in teratomic p oten tial is short range (as is the case for He atoms), then the minimum image con v en tion can b e used, whereb y a cuto distance of half the b o x length is used in the force calculations. This means that for a giv en He atom, only the neigh b oring atoms within a half b o xlength con tribute to the force on that He atom. In this w a y ev ery atom in the cen tral b o x con tributes to the force on ev ery other atom precisely once. If t w o atoms are within a half b o xlength of one another, then their in teraction is computed b et w een the corresp onding real atoms. If they are more than a half b o xlength from one another, then their in teraction o ccurs b et w een the real atom and the corresp onding virtual atom. In our sim ulations, w e b egan with a cen tral square b o x of length L = 40 : 46 au. This giv es a densit y = 0 : 00326 au 3 corresp onding to the kno wn bulk liquid He densit y 118 13 As for the initial state, w e w ould lik e to sim ulate the bulk He liquid in the canonical ensem ble, k eeping the n um b er of atoms, v olume and temp erature constan t. W e can pro vide a reasonable starting conguration b y randomly p ositioning the atoms within the cen tral b o x so that they are some minimal distance from one another. W e can also imp ose an initial temp erature b y assigning their momen ta acccording to the Maxw ellian distribution, p / exp p 2 2 mk B T ; (6.1) where k B is the Boltzmann constan t. W e p opulated our cen tral b o x using the minimal distance R min = 5 : 6 au for the conguration and temp erature T = 0 : 5 K for the momen tum distribution.
PAGE 118
104 There are man y metho ds of classical propagation, although w e found it useful to emplo y the v elo cit y V erlet algorithm, whic h is accurate to O ( t 2 ) and is selfstarting. In our case, it has the additional adv an tage of corresp onding directly to the propagation sc heme used for the quasiclassical v ariables in the relaxanddriv e pro cedure. While the c hoice of initial conditions ma y b e reasonable, it is necessary to equilibrate the system b y ev olving the atoms, and p erio dically rescaling the momen ta to obtain the desired temp erature, p p r T desir ed T actual : (6.2) There are other approac hes to temp erature calibration, but this rescaling sc heme a v oids c hanging the tra jectory directions so that the dynamics are b etter preserv ed. 34 Another imp ortan t issue in the equilibration is the timestep. T o k eep n umerical stabilit y w e imp osed a maxim um distance d max b y whic h an y atom could mo v e in a giv en timestep. This w as done b y nding the particle with maxim um sp eed v max and setting the (dynamically c hanging) timestep, t = d max v max : (6.3) Of course, this ignores the second order con tribution to the step in the v elo cit y V erlet algorithm, whic h dep ends on the particle accelerations. Ho w ev er, b y k eeping d max sucien tly small w e can neglect these higher order con tributions. W e found d max = 0 : 01 au to deliv er stable results without unnecessarily prolonging the sim ulation. Finally w e need to determine the eectiv e HeHe in teratomic p oten tial to generate the correct particle distribution from classical sim ulations. A t higher temp eratures, the HeHe p oten tial from Aziz and co w ork ers, 110 V Az iz ( R ), pro vides a
PAGE 119
105 v ery accurate description of the HeHe in teraction. Its form w as describ ed in Chapter 5, and its parameters for helium are sho wn in T able 6{1 Ho w ev er, b ecause of the quan tal eects at lo w temp eratures, w e exp ect an eectiv e p oten tial V ef f ( R ) substan tially dieren t from V Az iz ( R ) to b e necessary for the classical sim ulation to describ e the bulk He distribution. This can b e seen in Figure 6{2 where w e plot the radial distribution function obtained using V Az iz in our classical sim ulation, and compare it to path in tegral results. 118 Not only is the true radial distribution substan tially rattened b y appro ximately an order of magnitude, but it is also shifted to the righ t b y roughly 1 au. W e found that b y adding a correction, V 0 ( R ), to the Aziz p oten tial, V ef f ( R ) = V Az iz ( R ) + V 0 ( R ) ; (6.4) w e w ere able to generate a similar radial distribution function to results from pathin tegral calculations. 118 The eectiv e p oten tial con v erges to the Aziz p oten tial at short distances, but has an attractiv e w ell whic h is shifted to the righ t and is signican tly more shallo w than V Az iz The correction term w as constructed b y taking the dierence b et w een a scaled, shifted Aziz p oten tial and the original Aziz p otential, and m ultiplying b y a sigmoidal function that is asymptotically zero at short distances, V 0 ( R ) = ( R )[ V Az iz ( R R M =R s ) V Az iz ( R )] ; (6.5) ( R ) = 1 1 + exp [ a ( R R )] : (6.6) The parameters for V 0 are sho wn in T able 6{2 The p oten tials are graphed in Figure 6{3 and the eectiv e p oten tial is sho wn with expanded axes in Figure 6{4 Note that while the eectiv e p oten tial is substan tially ratter than the original Aziz p oten tial, it nonetheless retains an attractiv e w ell. Sim ulations rev ealed that the radial distribution function is v ery sensitiv e to the heigh t and p osition of this w ell.
PAGE 120
106 T able 6{1. P arameters for the HeHe in teraction from Aziz ( V Az iz ). R m A D C 6 C 8 C 10 5.6125 0.34648 1.922 10.735 1.893 1.414 1.349 0.414 0.171 T able 6{2. P arameters for the correction to the HeHe in teraction ( V 0 ) a R R S 4.0 4.0 0.003 6.62 0 2 4 6 8 10 12 14 2 4 6 8 10 12 14 16 g(R)R (au) V Aziz classical V Aziz path integral V eff classical 0 0.5 1 1.5 2 4 6 8 10 12 14 Figure 6{2. Radial distribution functions for bulk liquid helium.
PAGE 121
107 0.0006 0.0004 0.0002 0 0.0002 0.0004 0.0006 0.0008 0.001 2 3 4 5 6 7 8 9 10 Energy (au)R (au) V Aziz V' V eff = V Aziz + V' Figure 6{3. Comparison of the Aziz p oten tial with the eectiv e form.
PAGE 122
108 2e07 1.5e07 1e07 5e08 0 5e08 1e07 1.5e07 2e07 4 5 6 7 8 9 10 Energy (au)R (au) V eff Figure 6{4. Eectiv e HeHe p oten tial.
PAGE 123
109 6.4.2 Liquid Helium Droplets Ha ving formed an acceptable mo del of liquid He, w e are no w in a p osition to generate a mo del of a He droplet. W e b egin with a w ellequilibrated p erio dic b o x of liquid He at = 0 : 00326 au 3 and T = 0 : 5 K, and extract a sphere of 100 He atoms. Although there are w eak v an der W aals forces holding these atoms together, one m ust remem b er that w e ha v e used an eectiv e p oten tial with a v ery shallo w w ell, so that ev en at these ultralo w temp eratures, the He atoms will ultimately escap e in to the gas phase. T o coun ter this, w e in tro duce the constraining p oten tial discussed in Chapter 5, whic h consists of a sigmoidal function near the edge of the cluster. The parameters of the sigmoidal function w ere adjusted un til w e obtained a stable cluster whose densit y distribution w as consisten t with path in tegral and diusion Mon te Carlo calculations. 118 13 F or the constraining function, V hol d ( R ) = a 1 + exp [ b ( R R hol d )] ; (6.7) w e w ere able to generate an acceptable pure He droplet using a = 5 10 4 b = 0 : 08 and R hol d = 80 : 0. In Figure 6{5 w e graph the constraining p oten tial. Figure 6{6 sho ws that the temp erature of the drop is prop erly calibrated to T = 0 : 5 K, and that ructuations are small. Finally Figure 6{7 displa ys the densit y prole from the cen terofmass of the cluster, comparing our results with path in tegral and diusion Mon te Carlo calculations.6.4.3 LithiumHelium In teractions The LiHe p oten tials are describ ed using the pseudop oten tial formalism discussed in the previous c hapter. The parameters are dened in T ables 6{3 to 6{5 F or the electronic basis, w e used the 5s5p4d basis set of Cartesian Gaussian atomic functions tak en from Czuc ha j and co w ork ers. 90 In the w ork b y Rey es, 109 adiabatic p oten tial curv es including the d symmetries w ere ev aluated as a function of basis
PAGE 124
110 0 5e05 0.0001 0.00015 0.0002 0.00025 0.0003 0.00035 0.0004 0.00045 0.0005 0 20 40 60 80 100 120 140 160 Energy (au)R (au) V hold Figure 6{5. Constraining p oten tial used to k eep He atoms from ev ap orating.
PAGE 125
111 0.485 0.49 0.495 0.5 0.505 0.51 0.515 0.52 0 2e+06 4e+06 6e+06 8e+06 1e+07 Temperature (K)time (au) Figure 6{6. T emp erature ructations of the He droplet o v er time.
PAGE 126
112 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0 5 10 15 20 25 30 35 40 r (R)R (au) V eff classical V Aziz PIMC V Aziz DMC bulk 4 He Figure 6{7. Helium densit y prole from the cen terofmass of the cluster.
PAGE 127
113 T able 6{3. P arameters for the e Li in teraction. c c l B l l 0.1915 0.831 0 5.786 1.276 1 1.065 1.607 T able 6{4. P arameters for the e He in teraction. d B q B 1 l i B l l 1.3843 2.3265 0.706 0.79 0 1 0.83 1.30 0 2 2.27 0.50 0 2 2.27 0.50 1 1 0.12 0.75 1 2 1.87 1.00 size, and con v ergence w as found with the 5s5p4d basis. This basis pro vides an excellen t description of the lo w er adiabatic energies for the LiHe pair, as sho wn in Figure 6{8 F or this reason, w e emplo y ed this basis set for the remainder of our LiHe cluster sim ulations. F rom Figure 6{8 w e see that the 2 s ground state is essen tially repulsiv e. In fact, it con tains a v ery shallo w w ell near 11 au (not sho wn), whic h is resp onsible for the lo ose binding of Li atoms to the surface of He clusters. Ho w ev er, its depth is only 2 K, so that it will easily ev ap orate from the cluster. W e also see that the rst t w o excited states of lithium, 2 p and 2 p are substan tially dieren t in c haracter. In fact, the 2 p is a degenerate state of p orbitals in the plane p erp endicular to the LiHe axis, while 2 p con tains the p orbital parallel to the LiHe axis. The electronic o v erlap in the 2 p results in the curv e sho wing repulsion at all distances, while the 2 p curv e has an attractiv e w ell with a minim um around 4 : 5 au. These dierences result in v ery dieren t LiHe dynamics for dieren t electronic states.
PAGE 128
114 T able 6{5. P arameters for the LiHe core in teraction. A b R D 27.26 2.29 1.93 0.2 0.15 0.1 0.05 0 3 4 5 6 7 8 9 10 11 12 Energy (au)R (au) 2s s 2p p 2p s 3s s Figure 6{8. Adiabatic energy for Li and He as a function of in tern uclear distance.
PAGE 129
115 W e w ere also in terested in the eect of additional He atoms on the adiabatic energy giv en our in ten t to in tro duce the Li atom in to a cluster of He atoms. T o explore this eect, w e placed a single He atom at the origin, and then observ ed the c hange in adiabatic curv es as w e placed additional He atoms in the vicinit y of the origin. F or eac h conguration of helium atoms, w e mo v ed the lithium atom to w ard the origin along the zaxis, plotting the adiabatic energy as a function of the distance from the lithium atom to the origin. Initially w e added a single He atom at (0 ; 0 ; 6) au, and then another at (0 ; 0 ; 12) au. W e used a spacing of 6 au, since the radial distribution function p eak ed around this distance. The resulting adiabatic curv es are sho wn in Figure 6{9 where it is clear that the addition of the nal He atom aects the energy curv es negligibly In a second trial, w e b egan with a single He atom at the origin, and then added t w o along the yaxis: (0 ; 6 ; 0) au and (0 ; 6 ; 0) au. W e then added t w o more along the same line, (0 ; 12 ; 0) au and (0 ; 12 ; 0) au. The adiabatic curv es are plotted in Figure 6{10 where w e see the eects of additional He atoms are greater, but nonetheless the addition of the second set con tributes negligibly Finally w e lo ok ed at appro ximations to a cluster surface b y studying a 3 3 2 lattice of He atoms spanning [ 6 ; 6] [ 6 ; 6] [ 6 ; 0] au 3 with spacing of 6 au along eac h axis. The resulting adiabatic curv e w as compared to a larger 5 5 3 lattice of He atoms spanning [ 12 ; 12] [ 12 ; 12] [ 12 ; 0] au 3 again with 6 au spacing. The curv es are compared to the LiHe pair in Figure 6{11 The results rerect our exp ectations from the earlier curv es, in that there is v ery little dierence b et w een the smaller and larger surface in terms of adiabatic energies. These results giv e us condence that one can safely use a cuto in the order of 24 au (or ev en less) when computing the LiHe in teractions. This is esp ecially imp ortan t in larger cluster sim ulations, where the dominan t computational time is sp en t on matrix calculations.
PAGE 130
116 0.2 0.15 0.1 0.05 0 3 4 5 6 7 8 9 10 11 12 Energy (au)R (au) 2s s 2p p 2p s One He Two He Three He 0.13176 0.13168 0.1316 4.2 4.3 4.4 4.5 4.6 4.7 4.8 2p p 2s s 2p p 2p s Figure 6{9. Adiabatic energies for Li and one or more He along the zaxis.
PAGE 131
117 0.2 0.15 0.1 0.05 0 3 4 5 6 7 8 9 10 11 12 Energy (au)R (au) 2s s 2p p 2p s One He Three He Five He 0.1322 0.13204 0.13188 0.13172 4.2 4.3 4.4 4.5 4.6 4.7 4.8 2p p 2s s 2p p 2p s Figure 6{10. Adiabatic energy for Li and one or more He along the yaxis.
PAGE 132
118 0.2 0.15 0.1 0.05 0 3 4 5 6 7 8 9 10 11 12 Energy (au)R (au) 2s s 2p p 2p s One He 3x3x2 He 5x5x3 He 0.198115 0.198095 0.198075 9.5 10 10.5 11 11.5 12 12.5 13 2s s 2s s 2p p 2p s Figure 6{11. Adiabatic energy for Li and a surface of He atoms parallel to the xy plane.
PAGE 133
119 6.5 Results: Lithium Inside the Helium Cluster Exp erimen tal and quan tum mec hanical structure calculations sho w that the Li atom in v ariably resides on the surface of the He cluster, rather than within its bulk. W e examined this phenomenon b y b eginning with a liquid helium droplet in thermal equilibrium at T = 0 : 5 K. Our droplet con tained 100 He atoms, and w as brough t to equilibrium through the v elo cit y rescaling describ ed earlier. A t this p oin t, the He atom closest to the cen terofmass of the cluster w as replaced with a Li atom in the molecular ground state. This Li atom w as tak en to ha v e no momen tum initially so that all subsequen t dynamics w ere a result of in teractions with the surrounding He atoms. W e sho w the ev olution of the LiHe 99 cluster in Figure 6{12 where snapshots of the xy plane are compared at initial and nal times. As exp ected, the Li atom mo v es to w ard the cluster surface. What is in teresting, ho w ev er, is the sp eed at whic h it is ejected from the cluster's in terior. The path is direct, and motion from the cen ter to the surface (corresp onding to the t w o snapshots) tak es only 10,000 au (241 fs). This is a far greater sp eed than w ould b e exp ected through thermal motion at 0 : 5 K, and indicates that not only is the Li atom ejected from the cluster, but is ejected violen tly The con v erse is seen in exp erimen ts where Li can only b e forced in to liquid He through laser ablation. This is explained b y the strong repulsion of the ground state, esp ecially at t ypical distances separating He atoms, and sho ws that while the ground Li atom migh t b e pic k ed up on the cluster surface, it is extremely unlik ely to p enetrate this surface in to the cluster in terior. T o v erify that the motion of the Li atom is drastically dieren t from the thermal motion of the He atoms, w e p erformed the same sim ulation, but without replacing the He atom b y the Li atom. W e then follo w the ev olution of the cen tral He atom as it slo wly w anders through the cluster, and compare the motion to the Li atom in Figure 6{13 There are t w o ma jor dierences in their motion, as seen from the
PAGE 134
120 40 30 20 10 0 10 20 30 40 40 30 20 10 0 10 20 30 40 Y (au)X (au) A Initial 40 30 20 10 0 10 20 30 40 40 30 20 10 0 10 20 30 40 Y (au)X (au) B Initial Final Figure 6{12. Ev olution of ground state Li em b edded in the cen ter of a He cluster. A) Initial time t = 0 au. B) Final time t = 10,000 au.
PAGE 135
121 gure. The rst is that the He atom mak es essen tially a random w alk through the cluster, and sta ys within a cen tral v olume of radius 11 au. The second is that the He atom mo v es m uc h more slo wly than the Li atom. In fact, on the time scale of the Li motion, the He atom app ears nearly stationary whic h is wh y w e ha v e had to scale the He time axis b y t w o orders of magnitude to sho w an y appreciable mo v emen t. In Figure 6{14 w e trac k the electronic p opulations of the Li as it emerges from the cluster. W e see that the p opulation remains almost exclusiv ely in the ground state, with a ground p opulation loss of less than 0 : 5% b y the time the Li atom reac hes the surface. This small loss is lik ely attributable to the higher collision energies a v ailable to the Li atom as it gains momen tum, so that close approac hes to He atoms resulting in p opulation transfer are p ossible. One should note that at times earlier than 2000 au, b efore the Li atom pic ks up m uc h momen tum, p opulation transfer is not visibly detectable on the scale of the graph. This is a recurren t theme that w e will see again as w e ev olv e excited Li atoms, in that the v ery lo w kinetic energies are insucien t to bring the Li and He atoms close enough where curv e crossing, and therefore state mixing, are appreciable. 6.6 Results: Lithium on the Helium Cluster Surface Ha ving established that the Li atom do es not rest within the cluster, w e wish to in v estigate its b eha vior on the cluster surface. Sp ecically w e are in terested in follo wing the dynamics once the Li has b een excited to either the 2 p or 2 p state. Ho w ev er, since our cluster has a rather diuse surface, w e instead mo del the surface as a square lattice of He atoms. Reviewing the adiabatic curv es computed earlier, w e nd a lattice of 5 5 3 He atoms to b e sucien tly large to mo del a droplet surface. Using a square lattice also p ermits us to b etter dene the inital electronic state of the Li atom (2 p or 2 p ), and ev aluate the resp onse of the He atoms to the presence and motion of the attac hed Li atom.
PAGE 136
122 5 10 15 20 25 30 0 2000 4000 6000 8000 10000 distance from c.o.m. (au)time (au) Li(2s s ) He Figure 6{13. Comparison of Li and He motion within a He cluster. The time scale has b een reduced b y a factor of 100 for the He curv e.
PAGE 137
123 0.98 1 1.02 1.04 1.06 1.08 1.1 0 2000 4000 6000 8000 10000 h (2s s )time (au) 0.997 0.9975 0.998 0.9985 0.999 0.9995 1 1.0005 1.001 0 2000 4000 6000 8000 10000 Figure 6{14. Electronic p opulation of Li as it emerges from the He cluster.
PAGE 138
124 F or this lattice, w e nd a minim um energy along the ground state adiabatic curv e at R = 10 : 7 au, and th us place the Li atom initially at (0 ; 0 ; 10 : 7) au. Similar to the last section, w e b egin with the Li atom motionless, so that all n uclear motion comes from the LiHe in teractions. In the next t w o sections, w e follo w the dynamics of the Li atom once it has b een electronically excited. 6.6.1 Dynamics of Li( 2 p ) W e assume the Li atom has b een excited to 2 p and follo w its dynamics from this initial state. Snapshots of its ev olution are displa y ed in Figure 6{15 where w e see the motion of the He atoms and Li atom along the yz plane. Since this is a repulsiv e state, the Li atom mo v es a w a y from the He surface. W e also see the He atoms rep el in w ard from the Li atom. The snapshots are tak en at the initial time t = 0 au, and at the nal time t = 33 ; 000 au. Similar desorption from the surface is seen in exp erimen ts where either Na or K is excited to the repulsiv e state, and then immediately lea v es the surface. 11 The excited state deca ys slo wly and in the time frame of our sim ulation the Li atom remains in the 2 p state throughout. Ho w ev er, as sho wn in Figure 6{16 w e notice that while the total 2 p p opulation remains unit y there is almost complete transfer b et w een the 2 p and 2 p p opulations at around t = 2000 au. This is explained b y the fact that the 2 p state b ecomes triply degenerate at large distances, so that the 2 p and 2 p p opulations mix completely Ho w ev er, their degeneracy means that v ery little ligh t emission can b e exp ected. On the other hand, w e can induce dip ole emission b y p erturbing the system with an electromagnetic (EM) eld, ^ H 0 = ^ H ^ D E 0 cos ( t + ) ; (6.8) where ^ D is the system dip ole. F or an isolated Li atom, w e nd an energy dierence E = 0 : 068 au b et w een the ground and rst excited state. An EM eld with
PAGE 139
125 20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20 25 30 Y (au)Z (au) A Li 20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20 25 30 Y (au)Z (au) B Li Figure 6{15. Ev olution of Li(2 p ) as it recedes from the He cluster surface. A) Initial time t = 0 au. B) Final time t = 33 ; 000 au.
PAGE 140
126 0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 5000 10000 15000 20000 25000 30000 35000 htime (au) 2p p 2p s 2p p + 2p s Figure 6{16. Mixing of the Li(2 p ) and Li(2 p ) states at distances where Li(2 p ) is triply degenerate.
PAGE 141
127 frequency = 0 : 068 au 1 will stim ulate absorption and emission of ligh t b y the gas phase Li, as the Li atom exp eriences Li(2 p 2 s ) (called here the D ) transitions. As the EM frequency shifts from the electronic transition energy ligh t absorption and emission are still presen t, but with less magnitude. F or Li(2 p ) near the He surface, w e in tro duced an EM eld resonan t to the D line, with amplitude E 0 = 1 : 0 10 5 In addition, w e added the phase = 3 = 2 so that the electric eld rises from zero amplitude, rather than b eginning abruptly Finally w e c hose = (1 ; 1 ; 1) = p 3 so that there are comp onen ts along eac h axis, and the p erturbation has an eect regardless of the orien tation of the Li atom. In Figure 6{17 w e see the eect of in tro ducing the eld at t = 30 ; 000 au. By t = 32 ; 000 au, appro ximately half the total 2 p p opulation has deca y ed to the ground state. A t this p oin t, the 2 p p opulation b egins to rise again, corresp onding to ligh t absorption. This pattern of electronic deca y and reco v ery results from our semiclassical treatmen t of ligh tmatter in teractions. The eect of the EM eld is to induce the emission of dip ole radiation, whic h is measured from the second time deriv ativ e of the dip ole. Figure 6{18 sho ws the resulting dip ole sp ectrum for t w o time p erio ds, one when the p erturbing eld is applied during the initial 3000 au, and the second when the eld is applied during the nal 3000 au. W e also sho w the gas phase Li(2 p 2 s ) emission line, whic h w e compute from our asymptotic LiHe energies as 14903 cm 1 (and is extremely close to the kno wn v alue, 14904 cm 1 ). W e see that the Li(2 p ) emission p eak is initially blue shifted b y ab out 50 cm 1 ; at nal times, the p eak is substan tially broadened. This broadening can b e explained b y the motion of the He atoms, whic h ha v e acquired additional thermal energy from the repulsiv e in teraction with the Li atom.
PAGE 142
128 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 29000 30000 31000 32000 33000 h2ptime (au) without field with field Figure 6{17. Electronic p opulation of Li with and without a p erturbing electromagnetic eld, resonan t to the D line.
PAGE 143
129 13000 13500 14000 14500 15000 15500 16000 16500 17000 Intensity (arbitrary units)n (cm 1 ) initial time final time gas phase Li(2p <2s) Figure 6{18. Dip ole emission sp ectra of Li(2 p ) as it recedes from the He cluster surface.
PAGE 144
130 6.6.2 Dynamics of Li( 2 p ) The ev olution of the Li(2 p ) atom presen ts substan tially more in teresting features than the Li(2 p ) atom, as it is dra wn in to w ard the cluster surface. As suc h, it is not a matter of simply lea ving the cluster, but rather in v olv es non trivial in teractions with the He atoms. Snapshots are sho wn in Figure 6{19 again displa y ed as cross sections along the yz plane. Because of the attractiv e nature of the 2 p surface, w e see the Li atom indeed mo v es in to w ard the cluster surface. The Li atom ev en tually em b eds itself in the surface, and b y t = 67 ; 000 au it is completely surrounded b y He atoms. The He densit y around the Li atom eviden tly increases, although the HeHe repulsion coun ters this tendency once a certain He densit y is reac hed. Exp erimen tal results on excited Na atoms sho w the formation of NaHe N excimers, with N t ypically 3 or 4. According to the exp erimen tal results, once the excimers form, they desorb from the surface in 70 to 700 ps. 121 Our dynamics indicate that LiHe N excimers form on a m uc h shorter timescale (1.66 ps), alb eit w e are considering a dieren t alk ali atom. Ho w ev er, w e w ere not able to repro duce subsequen t desorption in longer sim ulation runs, partly b ecause of n umerical instabilities whic h ev en tually accum ulate. In Figure 6{20 w e sho w the p opulation of the 2 p state as the Li atom en ters the surface. P opulation transfer is exceptionally small, sho wing less than 0 : 01% o v er t = 67 ; 000 au, and once again reinforcing the conclusion that at these ultralo w temp eratures, collisioninduced p opulation transfer is negligible. Ho w ev er, similar to the study of Li(2 p ), w e are able to prob e the electronic state b y in tro ducing an electromagnetic eld and observing the dip ole emission sp ectrum. In Figure 6{21 w e sho w the dip ole sp ectrum resulting from an EM eld resonan t to the D line, in tro duced during the rst 3000 au of the sim ulation. The p eak of this sp ectrum o ccurs at = 14895 cm 1 whic h is shifted to the red of the gas phase emission b y 8 cm 1 Laserinduced ruorescence exp erimen ts sho w that a
PAGE 145
131 20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20 25 30 Y (au)Z (au) A Li 20 15 10 5 0 5 10 15 20 20 15 10 5 0 5 10 15 20 25 30 Y (au)Z (au) B Li Figure 6{19. Snapshot of Li(2 p ) as it in teracts with the He cluster surface. A) Initial time t = 0 au. B) Final time t = 67 ; 000 au.
PAGE 146
132 0.999 0.9992 0.9994 0.9996 0.9998 1 1.0002 1.0004 0 10000 20000 30000 40000 50000 60000 70000 h2p ptime (au) Figure 6{20. Electronic p opulation of Li(2 p ) as it in teracts with the He cluster surface.
PAGE 147
133 Li(2 s ) attac hed to a He cluster surface emits at a frequency whic h is red shifted b y ab out 10 cm 1 13 agreeing with our results. A m uc h more substan tial shift, ho w ev er, is found when the same EM eld is in tro duced in the nal 3000 au of the sim ulation, while the Li atom is surrounded b y the He atoms. Sho wn in Figure 6{22 the dip ole sp ectrum is shifted to the red b y appro ximately 4000 cm 1 While this v ery large shift initially app ears susp ect, it is explained b y examining the adiabatic energy curv es for a Li atom surr ounde d b y He atoms, rather than approac hing a surface. In this calculation, w e ha v e k ept the Li atom at the origin, and surrounded it with a square lattice of 8 He atoms (one at eac h corner of the lattice). W e ha v e then computed the adiabatic energy as a function of the distance the He atoms are lo cated along eac h axis. The curv es for the ground and rst excited state are sho wn in Figure 6{23 W e nd that ev en substan tially b efore the minim um of the w ell is reac hed, the distance b et w een the t w o curv es corresonds to the energy p eak of the dip ole emission sp ectrum. It is also in teresting to note ho w close the t w o curv es come to one another at small distances. This ma y explain wh y radiativ e deca y is not observ ed for Li in bulk He, where the n uclear congurations are forcibly suc h that the Li atom is alw a ys surrounded b y He atoms. In this case, it ma y b e energetically p ossible for excited Li atoms to reac h the left side of the adiabatic curv e, and deca y through curv e mixing in to the ground state. If this transfer w ere to o ccur b efore the time scale of sp on taneous emission, there w ould b e v ery little ligh t emission, and certainly none in the frequency range of the D emission line. In order to v alidate these sp ectral and adiabatic curv e results for Li(2 p ), w e in tro duced an electromagnetic eld resonan t with the supp osed sp ectral p eak, = 10 ; 730 cm 1 W e k ept the amplitude and phase the same as b efore, and turned on the eld at t = 64 ; 000 au. The results are sho wn in Figure 6{24 One can see that the 2 p p opulation deca ys steadily to ab out 25%, while the 2 s p opulation rises
PAGE 148
134 13000 13500 14000 14500 15000 15500 16000 16500 17000 Intensity (arbitrary units)n (cm 1 ) initial time gas phase Li(2p <2s) Figure 6{21. Dip ole emission sp ectrum of Li(2 p ) during the rst 3000 au.
PAGE 149
135 7000 8000 9000 10000 11000 12000 13000 14000 Intensity (arbitrary units)n (cm 1 ) Figure 6{22. Dip ole emission sp ectrum of Li(2 p ) during the nal 3000 au.
PAGE 150
136 0.2 0.18 0.16 0.14 0.12 0.1 2 2.5 3 3.5 4 4.5 5 5.5 6 Energy (au)R (au) D E = 10730 cm 1 2s s 2p p Figure 6{23. Adiabatic curv es of Li surrounded b y a cubic lattice of He atoms. The parameter R refers to the halflength of the lattice edge.
PAGE 151
137 corresp ondingly The EM eld frequency is rerected in the small ondulations, but the net p opulation transfer o v er time can b e in terpreted as induced ligh t emission. This sho ws that the classical electromagnetic eld can indeed b e used to induce b oth absorption and emission of ligh t. 6.7 Conclusion W e found that w e could repro duce the distribution of bulk liquid He at ultralo w temp eratures through a classical sim ulation b y using an eectiv e p oten tial whic h w as signican tly scaled do wn as w ell as shifted from the p oten tial prop osed b y Aziz and co w ork ers. Sp ecically w e w ere able to qualitativ ely repro duce the radial distribution function in the bulk. F urthermore, b y extracting a sphere of He atoms from this bulk, and imp osing a sigmoidal con tainmen t p oten tial on the system, w e w ere able to pro duce a liquid helium droplet whose densit y prole matc hed path in tegral and diusion Mon te Carlo results. Ha ving sucien tly mo deled the He droplet, w e in tro duced a ground state Li atom in to its cen ter, and follo w ed its dynamics. W e found that it w as rapidly exp elled from the cluster, and con trasted this motion from the slo w random w alk tak en b y a t ypical He atom. W e also found that the p opulation of the Li atom remained predominan tly in the ground state. W e then examined the dynamics of an excited Li(2 p ) atom near a mo del He surface. W e found that it mo v ed a w a y from the surface, while the He atoms recoiled somewhat in to the cluster. As with the Li(2 s ), the p opulation remained primarily in the original 2 p state. In order to prob e the electronic structure, w e in tro duced an electromagnetic eld to the system, resonan t to the gas phase Li(2 p 2 s ) line. W e computed the resulting dip ole sp ectrum, and found it to b e initially shifted appro ximately 50 cm 1 to the blue of the D line. A t later times, the p eak broadened, o wing to the increased thermal motion of the He atoms.
PAGE 152
138 0 0.2 0.4 0.6 0.8 1 64000 64500 65000 65500 66000 66500 67000 htime (au) Li(2p p ) Li(2s s ) Figure 6{24. Deca y of Li(2 p ) surrounded b y surface He atoms, induced b y an EM eld with frequency resonan t to the Li(2 p 2 s ) transition.
PAGE 153
139 W e also studied the dynamics of the excited Li(2 p ) atom near the mo del surface. In con trast to the 2 p state, the Li(2 p ) atom mo v ed to w ard the surface, and ev en tually found itself surrounded b y a high densit y of He atoms. These excimers ha v e b een sho wn to form for Na atoms near He clusters, follo w ed b y desorption of the excimer from the cluster. W e w ere able to sho w the formation of the excimers, but could not repro duce the subsequen t desorption, p ossibly due to the v ery long times required for the desorption to tak e place. An external electromagnetic eld ga v e rise to a dip ole emission sp ectrum initially shifted to the blue of the gas phase transition for Li(2 p ), and shifted to the red for Li(2 p ). These shifts are in agreemen t with exp erimen tal and other computational results. When the Li w as surrounded b y surface He atoms, the dip ole sp ectrum w as found to b e shifted to the red b y ab out 4000 cm 1 This v ery large shift w as explained b y lo oking at the adiabatic curv es of Li surrounded b y a matrix of He atoms, whic h rev ealed that the ground and rst excited states came sucien tly close to one another to corresp ond to the dip ole sp ectrum p eak. The adiabatic curv es also rev ealed ho w excited Li atoms in bulk He migh t b e able to deca y to the ground state with minimal ligh t emission, giv en the pro ximit y of the curv es. Finally w e in tro duced an electromagnetic eld resonan t with the dip ole emission p eak for Li(2 p ) surrounded b y He atoms. This electromagnetic eld induced deca y in to the ground state, un til the 2 s p opulation reac hed appro ximately 75%. This v alidated the sp ectral results, and demonstrated that the classical electromagnetic eld can b e used in mixed quan tumclassical sim ulations to induce b oth electronic excitation and deexcitation.
PAGE 154
CHAPTER 7 CONCLUSION The purp ose of this study has b een to dev elop a computationally feasible time dep enden t formalism to describ e the dynamics of fewand man yatom systems. W e ha v e dev elop ed a mixed quan tumclassical approac h whic h sho ws excellen t agreemen t with full quan tal calculations on small mo del systems, and is also capable of treating the ev olution of a realistic threedimensional mo del of an alk ali atom em b edded within a cluster of rare gas atoms. The general pro cedure used for the tests mo dels in tro duces an eectiv e p oten tial to the mixed quan tumclassical formalism, and uses classical tra jectories coupled to quan tal ev olution to generate the dynamics. The eciency of the calculations w as impro v ed b y in tro ducing a propagation algorithm whic h tak es in to accoun t the dieren t time scales of the n uclear and electronic motion, and pro vides a means to ev olv e the system with v arying timesteps, maximizing eciency while ensuring a prescrib ed lev el of accuracy In order to describ e the alk alirare gas clusters, the general formalism w as extended to incorp orate pseudop oten tial in teractions b et w een the alk ali v alence electron and the n uclear cores, and n umerical in terp olation metho ds w ere used to render the calculation of the Hamiltonian matrix elemen ts computationally feasible. In this conclusion, w e review the principal results of eac h c hapter. 7.1 Eectiv e P oten tial Quan tumClassical Liouville Equation Beginning with the quan tum Liouville equation of motion for the densit y op erator, w e divided the (general) system in to quan tal and quasiclassical degrees of freedom. By taking the partial Wigner transform o v er the quasiclassical degrees of freedom only and appro ximating for a large quasiclassical/quan tal mass ratio, 140
PAGE 155
141 w e deriv ed a mixed quan tumclassical equation of motion for the partially Wigner transformed densit y op erator. The transformed densit y op erator b ecomes a function of the quasiclassical phase space, but remains an op erator in quan tal space. By in tro ducing an eectiv e p oten tial to the mixed quan tumclassical equations of motion, w e generate a simplied v ersion, the eectiv e p oten tial quan tumclassical Liouville equation (EPQCLE), whic h can b e solv ed b y follo wing tra jectories through the quasiclassical phase space. This sc heme is rendered n umerically accessible b y in tro ducing a general basis set for the quan tal space, and b y discretizing the quasiclassical space on to a grid. The grid p oin ts ev olv e through quasiclassical tra jectories whic h are coupled to the quan tal state through the eectiv e force. The quan tal state ev olv es sim ultaneously using the relaxanddriv e metho d, whic h divides the quan tal motion in to a relaxation term com bined with a driving con tribution from the ev olving quasiclassical v ariables. The relaxanddriv e sc heme w as mo died to pro vide a robust v ariable timestep whic h accoun ted for the t ypically fast quan tal motion compared to the slo w er quasiclassical tra jectory ev olution. In the EPQCLE sc heme, measured quan tities w ere computed through an in tegration o v er quasiclassical phase space of the trace o v er quan tal v ariables. In this w a y exp ectation v alues of a general op erator, as w ell as p opulations of quan tal states, could b e computed analogously to the full quan tal calculations. 7.2 OneDimensional Tw oState Mo dels By sim ulationg three simple t w ostate mo dels, w e w ere able to compare the EPQCLE metho d with the exact quan tum mec hanical solution, obtained through the splitop erator fast F ourier transform grid metho d. W e found the EPQCLE to sho w excellen t agreemen t with the full quan tal results for man y observ ables, and at least qualitativ e agreemen t for all mo dels and energies. W e follo w ed a mo del of an alk ali atom aproac hing a metal surface, and ev aluated probabilit y transitions, p osition exp ectation and its deviation, momen tum
PAGE 156
142 exp ectation and its deviation, the densit y function and phase space grid ev olution. W e found the EPQCLE metho d repro duced qualitativ ely the exact v alues obtained through the quan tal treatmen t for all observ ables, and sho w ed the phase space grid to deform due to the action of the eectiv e p oten tial. W e examined a collision mo del in v olving t w o a v oided crossings, and found the transmission probabilit y obtained b y the EPQCLE to deviate from the quan tal results at lo w er energies. This problem w as observ ed in other theoretical mo dels whic h did not incorp orate n uclear coherence, and w as exp ected to arise b ecause of the imp ortance of n uclear coherence in the dual crossing mo del. Ho w ev er, w e found excellen t transmission probabilit y agreemen t at higher energies. Finally w e sim ulated a mo del of so dium in teracting with io dine, where the NaI complex oscillated b et w een a long range attractiv e ionic state and an asymptotically free neutral state. Disso ciation w as seen to o ccur at the diabatic crossing, as some of the ionic state w ould cross in to the neutral state and propagate to w ard free Na and I atoms. The oscillatory b eha vior w as repro duced qualitativ ely b y the EPQCLE mo del, and coherence b et w een the states sho w ed agreemen t with the quan tal results for long sim ulation times. Ev aluation of the phase space grid rev ealed b oth the oscillatory motion of the b ound NaI complex, and the disso ciated neutral state of the free atoms. 7.3 Alk aliRare Gas Clusters W e ha v e applied the EPQCLE formalism to the study of alk ali atoms em b edded within rare gas clusters. This has b een done b y treating the AkRg in teraction as a threeb o dy problem comp osed of the alk ali core, its v alence electrons and the rare gas core. In the con text of the EPQCLE, w e ha v e treated the electron as a quan tal v ariable and the atomic cores as quasiclassical. The electroncore in teractions ha v e b een describ ed using semilo cal l dep enden t pseudop oten tials, while the AkRg core in teractions ha v e b een describ ed with a tted parametric curv e. The
PAGE 157
143 RgRg in teraction w as mo deled b y a parametric curv e, and the threeb o dy system w as generalized to the case of one Ak atom and sev eral Rg atoms. In order to propagate the system, w e in tro duced a basis of segmen ted contractions of Gaussian atomic functions cen tered on the alk ali core. W e ev aluated all required matrix elemen ts explicitly and describ ed in detail the n umerical ev aluation of the HellmannF eynman force guiding the atomic cores. In addition, w e ha v e implemen ted a sc heme in v olving lo okup tables, in terp olation and rotation to ecien tly compute the pseudop oten tial matrix elemen ts, a ma jor b ottlenec k in the n umerical propagation of the cluster. W e ha v e studied the dynamics and sp ectra of a lithium atom em b edded in a cluster of helium atoms at v ery lo w temp eratures ( 0 : 5 K). W e follo w ed the path of the ground state lithium atom, and found it to rapidly mo v e from the cen ter of the cluster to its exterior. This is in agreemen t with exp erimen tal results whic h indicate that alk ali atoms reside on the surface of rare gas clusters and not within their in terior. W e ha v e also examined the dynamics of excited lithium atoms near a mo del helium surface, and found that the b eha vior is strongly dep enden t on the symmetry of the excited state. In the case of Li(2 p ), the atom is immediately rep elled from the surface, while the He atoms recoil sligh tly On the other hand, Li(2 p ) is attracted to the surface, and ev en tually em b eds itself completely within the He atoms. F or b oth excited states, w e induce dip ole emission b y in tro ducing a p erturbing electromagnetic eld resonan t to the Li(2 p 2 s ) line. The computed dip ole sp ectra ha v e b een calculated for initial and nal times, where w e nd sp ectral shifts to generally matc h other theoretical and exp erimen tal results. In the case of Li(2 p ) surrounded b y surface He atoms, w e observ e a v ery high red shift in the sp ectrum, whic h can b e explained b y examining the pro ximit y of the ground and rst excited state adiabatic curv es for a lithium atom approac hed b y a surrounding lattice of
PAGE 158
144 He atoms. This curv e pro ximit y ma y also explain the absence of radiation detected from excited Li in bulk liquid He, with radiationless deca y b eing p ossible through state mixing at sucien tly lo w LiHe distances. Finally w e w ere able to induce deca y and reco v ery of Li(2 p ) b y applying an electromagnetic eld resonan t to the Li(2 p 2 s ) transition, conrming our sp ectral results as w ell as demonstrating the abilit y of a classical electromagnetic eld to induce these electronic state c hanges in the con text of a mixed quan tumclassical sim ulation. 7.4 Soft w are Dev elopmen t Ov er the course of this pro ject, w e ha v e dev elop ed a general program, cauldron from scratc h in F ortran 90, to sim ulate the dynamics of arbitrary systems whic h can b e cast in to the EPQCLE formalism. Cauldron mak es no assumptions ab out and puts no limitations on the system size, and has b een dev elop ed in a mo dular fashion with the system, propagation and prop ert y calculations designed orthogonally Indeed, all three simple onedimensional mo dels, as w ell as the full LiHe N cluster, w ere sim ulated with precisely the same relaxanddriv e propagation co de. This orthogonal approac h has enabled cauldron to b e extended as new questions arose and new systems w ere studied, without the need to mo dify already w orking co de. W e ha v e also dev elop ed a second program, qualdron analgous to cauldron but for sim ulating onedimensional t w ostate systems using the splitop erator fast F ourier transform formalism. Qualdron has b een used primarily to generate precise results for comparison with the ( cauldron ) EPQCLE results. Qualdron has also b een dev elop ed with orthogonal mo dules for the system, propagation and prop erties, but is limited to onedimensional t w ostate systems describ ed in a diabatic basis. Ho w ev er, the extension to m ultistate systems w ould b e a trivial task, requiring the mo dication of a couple of subroutines in the propagation mo dule. Both cauldron and qualdron are large programs, in v olving o v er 35 ; 000 lines of co de and 200 subroutines. Both their use and in ternal functioning ha v e b een
PAGE 159
145 extensiv ely describ ed in t w o accompan ying man uals, where sp ecic proto cols for extending the co des to new problems ha v e b een suggested. 7.5 F uture W ork While w e ha v e fo cused on the application of the EPQCLE to three small mo dels as w ell as the realistic threedimensional LiHe N cluster, it w ould b e in teresting to ev aluate the p erformance and accuracy of the EPQCLE for additional test mo dels where precise quan tal results can b e computed, in order to further assess the b enets and limitations of the metho d. It w ould also b e in teresting to examine the dynamics of AkRg N clusters in v olving hea vier alk ali or rare gas atoms. Finally it w ould b e instructiv e to study the dynamics and sp ectra of collisions of Ak atoms with Rg atoms at higher temp eratures, where the Rg clusters v ap orize in to the gas phase. AkRg collisions at thermal and h yp erthermal energies ha v e b een studied recen tly b y Rey es, and cauldron pro vides a means to study collisions at these energies where m ultiple rare gas atoms are in v olv ed.
PAGE 160
APPENDIX A THE CAULDRON PR OGRAM A.1 Ov erview Cauldron has b een written to n umerically solv e the EPQCLE for an arbitrary n um b er of quan tal and classical degrees of freedom. F rom the programmer's p oin t of view, cauldron solv es the follo wing partial dieren tial matrix equation: @ @ t = a + b + ~ c @ @ ~ R + ~ d @ @ ~ P ; (A.1) where the matrices a b and ha v e dimensions of the electronic basis set, and the v ectors ~ c ~ d ~ R and ~ P ha v e dimensions rerecting the classical degrees of freedom. is a function of R P and t while the co ecien t v ectors and matrices can b e functions of these v ariables as w ell as itself (for example, when the HellmannF eynman eectiv e force is used). P osed this w a y Eq. A.1 can represen t man y systems and can b e solv ed using a v ariet y of propagation metho ds. T o this end, cauldron w as designed along three orthogonal branc hes, reminiscen t of classes used in ob jectorien ted design. Ho wev er, while v ery mo dular, cauldron is not truly ob jectorien ted in either design or implemen tation, and in the in terest of a v oiding confusion and misrepresen tation, ob jectorien ted terminology will not b e used. The three orthogonal branc hes of cauldron are: 1. System: The system branc h pro vides the time dep enden t co ecien ts for Eq. A.1 completely dening the mo del under study All Hamiltonian terms (quan tal, classical and coupling), as w ell as the eectiv e p oten tial, are con tained within the co ecien ts and ev aluated en tirely within the system branc h. New mo dels 146
PAGE 161
147 are implemen ted b y dev eloping a system branc h sp ecic to the mo del, without reference to an y other branc h. 2. Propagation: The propagation branc h uses the co ecien ts computed b y the system branc h to ev olv e the system. This branc h is general, and a giv en metho d (e.g., relaxanddriv e) ma y b e applied to all systems without mo dication, regardless of system dimension. 3. Prop erties: A t eac h timestep, the prop erties branc h p erforms calculations on the system (e.g. exp ectation v alues of observ ables) and sends the results to appropriate output les for p ostpro cessing. The main cauldron executable calls the appropriate co de from eac h branc h, as sho wn in Figure A{1 Since these co de branc hes are orthogonal to one another, it is a simple matter to add a new system, propagation metho d or prop ert y calculation of arbitrary complexit y without aecting the rest of the co de. A.2 Comp onen t Descriptions A.2.1 Read Input File Eac h sim ulation is describ ed b y a single input le, cauldron.init comp osed of three namelists of data: one for the system, one for the propagation, and one for the prop ert y calculations. T o run a sim ulation, cauldron m ust b e executed in a directory con taining the desired cauldron.init le. A.2.2 System: Get Dieren tial Equation Co ecien ts A t this p oin t, the system branc h pro vides the dieren tial equation co ecien ts at a giv en time. Curren tly implemen ted systems are: Alk ali atommetal surface. This is a onedimensional t w ostate mo del of an alk ali atom approac hing a metal surface. Dual a v oided crossing. This system mo dels a collision b et w een t w o n uclei in v olving t w o electronic states and t w o a v oided crossings.
PAGE 162
148 Evolve single timestep Propagation: yes no end begin finished? Properties: Output properties Get diffeq coefficients System: read input file Figure A{1. Flo w c hart describing the cauldron program.
PAGE 163
149 NaI complex. Here w e mo del the NaI complex with the n uclear separation and t w o diabatic electronic surfaces. Bulk liquid helium. This is a threedimensional mo del of an arbitrary n um b er of classical helium atoms in a p erio dically replicating cen tral v olume. Helium droplet. Here w e mo del a liquid droplet of classically in teracting helium atoms. Lithiumhelium cluster. This system mo dels the in teraction of a lithium atom within a cluster of helium atoms in the framew ork of the EPQCLE, treating the atomic cores quasiclassically and the lithium v alence electron quan tally Their in teraction is describ ed through l dep enden t semilo cal pseudop oten tials and the eectiv e HellmannF eynman force. A.2.3 Propagation: Ev olv e Single Timestep The propagation branc h uses the the dieren tial equation co ecien ts to adv ance the system a single timestep. Curren tly implemen ted branc hes are: Relaxanddriv e. This is the mixed quan tumclassical relaxanddriv e algorithm, with the capabilit y to use either a xed or v ariable timestep. The propagation can also b e split b et w een the relaxation and the driving terms, omitting the driving term if desired. V elo cit y V erlet with p erio dic b oundaries. This is a classical v elo cit y V erlet propagation, but using p erio dic b oundary conditions so that particles lea ving the cen tral v olume reapp ear on the opp osite side. V elo cit y V erlet without p erio dic b oundaries. This is the classical v elo cit y V erlet algorithm with innite b oundaries, useful for propagating isolated clusters of helium. A.2.4 Prop erties: Output Prop erties The prop erties branc h outputs v arious quan tities of in terest to sp ecic and predened output les. Curren tly implemen ted prop erties are:
PAGE 164
150 Exp ectation v alue of p osition and momen tum. P osition and momen tum disp ersion. T race (b oth quan tum and classical) of the partially Wigner transformed densit y matrix. State probabilities. Nuclear conguration, for b oth original and smo othed surfaces. T emp erature, kinetic and p oten tial energy Nuclear densit y prole. P air distribution function. Dip ole exp ectation v alues. A.3 Subroutine Details Cauldron has b een written en tirely from scratc h in F ortran 90, and in v olv es nearly 10,000 lines of co de in o v er 100 les. It w ould b e a hop eless task to describ e the details in the space of an app endix. Ho w ev er, a detailed reference man ual 123 has b een written to accompan y cauldron whic h guides the programmer through all asp ects of b oth use and extension of the co de.
PAGE 165
APPENDIX B SPLIT OPERA TORF AST F OURIER TRANSF ORM METHOD Consider a quan tal system describ ed b y ( r ; R ; t ). The BornOpp enheimer separation builds this w a v efunction in terms of a complete and orthonormal set f k ( r ; R ) g of electronic functions, ( r ; R ; t ) = X l l ( r ; R ) ( R ; t ) : (B.1) Beginning with the time dep enden t Sc hr odinger equation, i ~ @ @ t = ^ H ; (B.2) w e can insert a general Hamiltonian in v olving n uclear (N) and electronic (e) co ordinates, ^ H = ^ K N + ^ T e + ^ V N N + ^ V N e + ^ V ee  {z } ^ H e (B.3) so that i ~ @ @ t = [ ^ K N + ^ H e ] : (B.4) Switc hing to Dirac notation for con v enience, and setting ~ = 1, w e m ultiply on the left b y h k j h k j i @ @ t X l j l ij l i = h k j ^ K N + ^ H e ( j l ij l i ) (B.5) so that i @ @ t j k i = h k j ^ K N ( j l ij l i ) + h k j ^ H e j l i  {z } ^ V k l j l i : (B.6) 151
PAGE 166
152 Expanding the rst term on the RHS, h k j ^ K N ( j l ij l i ) = h k j 1 2 M r 2R ( j l ij l i ) = h k j 1 2 M r 2R j l ij l i + j l ir 2R j l i : (B.7) In the diabatic basis, the momen tum coupling h k jr R j l i disapp ears, so that i @ @ t k = 1 2 M r 2R k + ^ V k l l : (B.8) F or a t w ostate system, denoting the column v ector ( 1 ; 2 ), w e can write @ @ t = i 0B@ 1 2 M r 2R 0 0 1 2 M r 2R 1CA  {z } ^ K i 0B@ ^ V 11 ^ V 12 ^ V 21 ^ V 22 1CA  {z } ^ V : (B.9) F ormally w e can ev olv e a timestep t using the timeordered op erator ^ T : ( t ) = ^ T exp 24 t Z 0 ( i ^ K ( t 0 ) i ^ V ( t 0 )) dt 0 35 (0) : (B.10) If the v ariation of the Hamiltonian is small on the timescale of t w e can write Eq. B.10 as ( t ) = exp h ( i ^ K (0) i ^ V (0)) t i (0) : (B.11) The split op erator metho d breaks up the exp onen tial, ( t ) = exp i t 2 ^ K (0) exp h i t ^ V (0) i exp i t 2 ^ K (0) (0) + O ( t 3 ) : (B.12)
PAGE 167
153 Since ^ K is diagonal but nonlo cal in the n uclear co ordinate space, w e can ev alute its exp onen tial b y computing the deriv ativ es in F ourier space, ~ ( k ) = r 1 2 1 Z 1 ( R ) e ik R dR ) @ 2 @ R 2 ~ ( k ) = k 2 ~ ( k ) On the other hand, ^ V is lo cal but nondiagonal. T o ev aluate its exp onen tial, w e need to diagonalize ^ V : ^ V = D 1 D ) exp [ i t ^ V ] = D exp[ i t ] D 1 = D 0B@ exp ( i t 1 ) 0 0 exp ( i t 2 ) 1CA D 1 : Th us, the split op erator fast F ourier transform metho d ev olv es the quan tal system a single timestep as follo ws: 1. Ev aluate exp [ i t= 2 ^ K ] ( t ) b y computing the second spatial deriv ativ e in F ourier space. 2. Apply exp [ i t ^ V ( t )] b y diagonalizing ^ V ( t ) and then con v erting the exp onential of the matrix in to a matrix of exp onen tials. 3. Apply exp [ i t= 2 ^ K ] once again using the fast F ourier transform. In practice, when sev eral timesteps are computed in a ro w, adjoining exp onen tials of ^ K (i.e., step 1/ follo wing step 3/) can b e com bined in to a single op eration with double the timestep.
PAGE 168
APPENDIX C THE QUALDRON PR OGRAM C.1 Ov erview Qualdron has b een written to n umerically solv e the time dep enden t Sc hr odinger equation for a t w ostate electronic system with a single classical degree of freedom. F rom the programmer's p oin t of view, qualdron solv es the follo wing partial dieren tial equation: i @ @ t = ( ^ K + ^ V ) ; (C.1) where is the column v ector of t w o electronic states, = ( 1 ; 2 ). The basis is assumed to b e diabatic, so that ^ K is diagonal and ^ V con tains the coupling b et w een states. While qualdron has b een written for the t w ostate case, expansion to an arbitrary n um b er of states w ould b e straigh tforw ard, and require the replacemen t of a single subroutine to rerect the m ultidimensional problem. P osed this w a y Eq. C.1 can represen t man y dieren t systems b y c hanging the v alue of the p oten tial coupling elemen ts ^ V in the Hamiltonian. T o this end, qualdron w as designed along three orthogonal branc hes, similar to the design of cauldron These branc hes are: 1. System: The system branc h pro vides the initial v alues and Hamiltonian for Eq. C.1 This section of co de denes the system, and can represen t an y onedimensional t w ostate mo del that can b e describ ed b y a Hamiltonian in the diabatic represen tation. 2. Propagation: The propagation branc h uses the Hamiltonian computed b y the system branc h to ev olv e the system. The propagation co de is general, and 154
PAGE 169
155 a giv en metho d (e.g., split op erator) ma y b e applied to all dened systems without mo dication. 3. Prop erties: A t eac h timestep, the propagation co de calls the prop erties co de, whic h mak es system calculations and outputs appropriate les. The main qualdron executable calls the appropriate co de from eac h area to create the sim ulation, as sho wn in Figure C{1 Since these co de branc hes are for the most part orthogonal to one another, it is a simple matter to add a new system, a new propagation metho d, or new prop ert y calculations without aecting the rest of the co de. C.2 Comp onen t Descriptions C.2.1 Read Input File Eac h sim ulation is describ ed b y a single input le, qualdron.init comp osed of three namelists of data: one for the system, one for the propagator, and one for the prop erties. T o run a sim ulation, qualdron m ust b e executed in a directory con taining the desired qualdron.init le. C.2.2 System: Get Hamiltonian Matrix Elemen ts A system pro vides the Hamiltonian matrix elemen ts at a giv en time. Curren tly implemen ted systems are: Alk ali atommetal surface. Dual a v oided crossing. NaI complex. C.2.3 Propagation: Ev olv e Single Timestep The propagation co de uses the Hamiltonian elemen ts to adv ance the system a single timestep. Curren tly the only implemen ted propagation is the split op erator fast F ourier transform metho d.
PAGE 170
156 Propagation: Evolve single trajectory Get Hamiltonian System: yes no end begin finished? Properties: Output properties read input file Figure C{1. Flo w c hart describing the qualdron program.
PAGE 171
157 C.2.4 Prop erties: Output Prop erties The prop erties co de outputs v arious quan tities of in terest to sp ecic and predened output les. Curren tly implemen ted prop erties are: Exp ectation v alue of p osition and momen tum. P osition and momen tum disp ersion. State probabilities. Nuclear congurations. C.3 Subroutine Details Qualdron has b een written en tirely from scratc h in F ortran 90, and in v olv es nearly 2,500 lines of co de in o v er 100 les. Since these details could not b e adequately describ ed b y an app endix, a detailed reference man ual 124 has b een written to accompan y qualdron whic h guides the programmer through all asp ects of b oth use and extension of the co de. In particular, instructions for c hanging qualdron from a t w ostate to a m ultistate system are presen ted in detail.
PAGE 172
REFERENCE LIST [1] Aatto Laaksonen and Y ao quan T u. Mole cular Dynamics: F r om Classic al to Quantum Metho ds v olume 7 of The or etic al and Computational Chemistry c hapter 1. Elsevier, Amsterdam, 1999. [2] J. Bric kmann and U. Sc hmitt. Mole cular Dynamics: F r om Classic al to Quantum Metho ds v olume 7 of The or etic al and Computational Chemistry c hapter 2. Elsevier, Amsterdam, 1999. [3] Nancy Makri. Timedep enden t quan tum metho ds for large systems. A nnual R eview of Physic al Chemistry 50:167{191, 1999. [4] M. Hillery R. F. O'Connell, M. O. Scully and E. P Wigner. Distribution functions in ph ysics: F undamen tals. Physics R ep orts 106(3):121{167, 1984. [5] Chiac hin Tso o, Dario A. Estrin, and Sherwin J. Singer. Electronic energy shifts of a so dium atom in argon clusters b y sim ulated annealing. Journal of Chemic al Physics 93(10):7187{7200, 1990. [6] Kenneth Haug and Horia Metiu. The absorption sp ectrum of a p otassium atom in a Xe cluster. Journal of Chemic al Physics 95(8):5670{5680, 1991. [7] Kenneth Haug and Horia Metiu. Absorption sp ectrum calculations using mixed quan tumGaussian w a v e pac k et dynamics. Journal of Chemic al Physics 99(9):6253{6263, 1993. [8] Glenn Mart yna, Ching Cheng, and Mic hael L. Klein. Electronic states and dynamical b eha vior of LiXe n and CsXe n clusters. Journal of Chemic al Physics 95(2):1318{1336, 1991. [9] F. Stienk emeier, J. Higgins, C. Callegari, S. I. Kanorsky W. E. Ernst, and G. Scoles. Sp ectroscop y of alk ali atoms (Li, Na, K) attac hed to large helium clusters. Zeitschrift f ur Physik D 38:253{263, 1996. [10] Y ongkyung Kw on and K. Birgitta Whaley A tomicscale quan tum solv ation structure in sup erruid helium4 clusters. Physic al R eview L etters 83(20):4108{ 4111, 1999. [11] F rank Stienk emeier and Andrey F. Vileso v. Electronic sp ectroscop y in He droplets. Journal of Chemic al Physics 115(22):10119{10137, 2001. [12] R. N. Barnett and K. B. Whaley Molecules in helium clusters: SF 6 He N Journal of Chemic al Physics 99(12):9730{9744, 1993. 158
PAGE 173
159 [13] Akira Nak a y ama and Noic hi Y amashita. P ath in tegral Mon te Carlo study on the structure and absorption sp ecta of alk ali atoms (Li, Na, K) attac hed to sup erruid helium clusters. Journal of Chemic al Physics 114(2):780{790, 2001. [14] J. P eter T o ennies and Andrei F. Vileso v. Sp ectroscop y of atoms and molecules in liquid helium. A nnual R eview in Physic al Chemistry 49:1{41, 1998. [15] P W. A tkins and R. S. F riedman. Mole cular Quantum Me chanics Oxford Univ ersit y Press, Oxford, 3rd edition, 1997. [16] Ro dney J. Bartlett, editor. R e c ent A dvanc es in Couple d Cluster Metho ds v olume 3 of R e c ent A dvanc es in Computational Chemistry W orld Scien tic, Singap ore, 1997. [17] Mic hael Springb org. Metho ds of Ele ctr onicStructur e Calculations: F r om Mole cules to Solids Wiley New Y ork, 2000. [18] A. Szab o and N. S. Ostlund. Mo dern Quantum Chemistry Do v er, Mineola, New Y ork, 1996. [19] Douglas L. Strout and Gusta v o E. Scuseria. A quan titativ e study of the scaling prop erties of the HartreeFo c k metho d. Journal of Chemic al Physics 102(21):8448{8452, 1995. [20] Da vid J. Griths. Intr o duction to Quantum Me chanics Pren tice Hall, Upp er Saddle Riv er, New Jersey 1994. [21] Karl Blum. Density Matrix The ory and Applic ations Plen um Press, New Y ork, 2nd edition, 1981. [22] John D. Simon, editor. Ultr afast Dynamics of Chemic al Systems v olume 7 of Understanding Chemic al R e activity Klu w er Academic Publishers, Dordrec h t, 1994. [23] R. Guan tes and S. C. F aran tos. High order nite dierence algorithms for solving the Sc hro dinger equation in molecular dynamics. Journal of Chemic al Physics 111(24):10827{10835, 1999. [24] D. A. Maziotti. Sp ectral dierence metho ds for solving dieren tial equations. Chemic al Physics L etters 299(5):473{480, 1999. [25] Mark Thac h uk and George C. Sc hatz. Timedep enden t metho ds for calculating thermal rate co ecien ts using rux correlation functions. Journal of Chemic al Physics 97(10):7297{7313, 1992. [26] P atric k L. Nash and L. Y. Chen. Ecien t nite dierence solutions to the timedep enden t Sc hr odinger equation. Journal of Computational Physics 130:266{ 268, 1997.
PAGE 174
160 [27] D. Koslo and R. Koslo. A Fourier metho d solution for the time dep enden t Sc hro dinger equation as a to ol in molecular dynamics. Journal of Computational Physics 52(1):35{53, 1983. [28] T. W. K orner. F ourier A nalysis Univ ersit y Press, Cam bridge, 1988. [29] C. Zenger. Par al lel A lgorithms for Partial Dier ential Equations v olume 31 of Notes on Numeric al Fluid Me chanics c hapter Sparse Grids. View eg, Braunsc h w eig, 1991. [30] Dmitrii V. Shalashilin and Mark S. Child. Time dep enden t quan tum propagation in phase space. Journal of Chemic al Physics 113(22):10028{10036, 2000. [31] Gert D. Billing. Timedep enden t quan tum dynamics in a GaussHermite basis. Journal of Chemic al Physics 110(12):5526{5537, 1999. [32] Herb ert Goldstein. Classic al Me chanics AddisonW esley Massac h usetts, 2nd edition, 1980. [33] Harv ey Gould and Jan T ob o c hnik. A n Intr o duction to Computer Simulation Metho ds AddisonW esley Reading, Massac h usetts, 2nd edition, 1996. [34] M. P Allen and D. J. Tildesley Computer Simulation of Liquids Oxford Univ ersit y Press, Oxford, 1987. [35] J. Barnes and P Hut. A hierarc hical O(N log N) force calculation algorithm. Natur e 324:446{449, 1986. [36] L. Greengard and V. Rokhlin. A fast algorithm for particle sim ulations. Journal of Computational Physics 73:325{348, 1987. [37] Leslie Greengard. F ast algorithms for classical ph ysics. Scienc e 265:909{914, 1994. [38] Celeste Sagui and Thomas A. Darden. Molecular dynamics sim ulations in biomolecules: Longrange electrostatic eects. A nnual R eview of Biophysics and Biomole cular Structur e 28:155{179, 1999. [39] Da vid A. Mic ha. Timedep enden t man yelectron treatmen t of electronic energy and c harge transfer in atomic collisions. Journal of Physic al Chemistry A 103(38):7562{7574, 1999. [40] Da vid A. Mic ha. Timeev olution of m ulticonguration densit y functions driv en b y n uclear motions. International Journal of Quantum Chemistry 60:109{118, 1996. [41] Da vid A. Mic ha. T emp oral rearrangemen t of electronic densities in slo w atomic collisions. International Journal of Quantum Chemistry 51:499{518, 1994.
PAGE 175
161 [42] Lic hang W ang and Anne B. McCo y Timedep enden t studies of reaction dynamics: a test of mixed quan tum/classical timedep enden t selfconsisten t eld appro ximations. Physic al Chemistry and Chemic al Physics 1:1227{1235, 1999. [43] Eric J. Heller. Timedep enden t approac h to semiclassical dynamics. Journal of Chemic al Physics 62(4):1544{1555, 1975. [44] Eric J. Heller. Time dep enden t v ariational approac h to semiclassical dynamics. Journal of Chemic al Physics 64(1):63{73, 1976. [45] John C. T ully T ra jectory surface hopping approac h to nonadiabatic molecular collisions: the reaction of H + with D 2 Journal of Chemic al Physics 55(2):562{ 572, 1971. [46] E. Deumens, A. Diz, H. T a ylor, and Y. Ohrn. Timedep enden t dynamics of electrons and n uclei. Journal of Chemic al Physics 96(9):6820{6833, 1992. [47] E. Deumens, A. Diz, R. Longo, and Y. Ohrn. Timedep enden t theoretical treatmen ts of the dynamics of electrons and n uclei in molecular systems. 66(3):917{983, 1994. [48] Nancy Makri. F eynman path in tegration in quan tum dynamics. Computer Physics Communic ations 63:389{414, 1991. [49] R. F eynman and A. R. Hibbs. Quantum Me chanics and Path Inte gr als McGra wHill, New Y ork, 1965. [50] M. BenNun and T o dd J. Martinez. A m ultiple spa wning approac h to tunneling dynamics. Journal of Chemic al Physics 112(14):6113{6121, 2000. [51] John C. T ully Molecular dynamics with electronic transitions. Journal of Chemic al Physics 93(2):1061{1071, 1990. [52] John R. Klauder and BoSture Sk agerstam. Coher ent States: Applic ations in Physics and Mathematic al Physics W orld Scien tic, Singap ore, 1985. [53] La wrence S. Sc h ulman. T e chniques and Applic ations of Path Inte gr ation Wiley New Y ork, 1981. [54] Nancy Makri. Numerical path in tegral tec hniques for long time dynamics of quan tum dissipativ e systems. Journal of Mathematic al Physics 36(5):2430{ 2457, 1995. [55] Nancy Makri and Dmitrii E. Mak aro v. T ensor propagator for iterativ e quantum time ev olution of reduced densit y matrices. I I. n umerical metho dology Journal of Chemic al Physics 102(11):4611{4618, 1995.
PAGE 176
162 [56] Xiong Sun, Haobin W ang, and William H. Miller. Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized appro ximation to the initial v alue represen tation. Journal of Chemic al Physics 109(17):7064{ 7074, 1998. [57] Xiong Sun and William H. Miller. Mixed semiclassicalclassical approac hes to the dynamics of complex molecular systems. Journal of Chemic al Physics 106(3):916{927, 1997. [58] Xiong Sun and William H. Miller. Semiclassical initial v alue represen tation for electronically nonadiabatic molecular dynamics. Journal of Chemic al Physics 106(15):6346{6353, 1997. [59] U. F ano. Description of states in quan tum mec hanics b y densit y matrix and op erator tec hniques. R eviews of Mo dern Physics 29(1):74{93, 1957. [60] Bret Jac kson. Reduced densit y matrix description of gassolid in teractions: Scattering, trapping and desorption. Journal of Chemic al Physics 108(3):1131{1139, 1998. [61] Jonathan Grad, Yi Jing Y an, Azizul Haque, and Shaul Muk amel. Reduced equations of motion for semiclassical dynamics in phase space. Journal of Chemic al Physics 86(6):3441{3454, 1987. [62] Y ang Zhao, Satoshi Y ok o jima, and GuanHua Chen. Reduced densit y matrix and com bined dynamics of electrons and n uclei. Journal of Chemic al Physics 113(10):4016{4027, 2000. [63] Herman J. C. Berendsen and Janez Ma vri. Quan tum sim ulation of reaction dynamics b y densit y matrix ev olution. Journal of Physic al Chemistry 97:13464{13468, 1993. [64] Herman J. C. Berendsen and Janez Marvi. Approac h to nonadiabatic transitions b y densit y matrix ev olution and molecular dynamics sim ulation. International Journal of Quantum Chemistry 57:975{983, 1996. [65] Marc F. Lensink, Janez Ma vri, and Herman J. C. Berendsen. Sim ultaneous in tegration of mixed quan tumclassical systems b y densit y matrix ev olution equations using in teraction represen tation and adaptiv e time step in tegrator. Journal of Computational Chemistry 17(11):1287{1295, 1996. [66] Da vid A. Mic ha and Keith Runge. Timedep enden t man yelectron approac h to slo w ionatom collisions: The coupling of electronic and n uclear motions. Physic al R eview A 50(1):322{336, 1994. [67] Keith Runge and Da vid A. Mic ha. Timedep enden t approac h to slo w ionatom collisions for systems with one activ e electron. Physic al R eview A 53(3):1388{ 1399, 1996.
PAGE 177
163 [68] Keith Runge and Da vid A. Mic ha. Timedep enden t man yelectron approac h to slo w ionatom collisions for systems with sev eral activ e electrons. Physic al R eview A 62(2):22703{1{22703{10, 2000. [69] Keith H. Hughes and Rob ert E. Wy att. T ra jectory approac h to dissipativ e quan tum phase dynamics: Application to barrier scattering. Journal of Chemic al Physics 120(9):4089{4097, 2004. [70] Irene Burghardt and G erard P arlan t. On the dynamics of coupled Bohmian and phasespace v ariables: A new h ybrid quan tumclassical approac h. Journal of Chemic al Physics 120(7):3055{3058, 2004. [71] I. Burghardt and L. S. Cederbaum. Hydro dynamic equations for mixed quan tum states. I. General form ulation. Journal of Chemic al Physics 115(22):10303{10311, 2001. [72] I. Burghardt and L. S. Cederbaum. Hydro dynamic equations for mixed quan tum states. I I. Coupled electronic states. Journal of Chemic al Physics 115(22):10312{10322, 2001. [73] Jerem y B. Maddo x and Eric R. Bittner. Quan tum relaxation dynamics using Bohmian tra jectories. Journal of Chemic al Physics 115(14):6309{6316, 2001. [74] Da vid A. Mic ha and Brian Thorndyk e. Dissipativ e dynamics in man yatom systems: A densit y matrix treatmen t. International Journal of Quantum Chemistry 2002. T o b e published. [75] Arnoldo Donoso and Craig C. Martens. Sim ulation of coheren t nonadiabatic dynamics using classical tra jectories. Journal of Physic al Chemistry A 102:4291{4300, 1998. [76] Arnaldo Donoso and Craig C. Martens. Semiclassical m ultistate Liouville dynamics in the adiabatic represen tation. Journal of Chemic al Physics 112(9):3980{3989, 2000. [77] Arnaldo Donoso, Daniela Kohen, and Craig C. Martens. Sim ulation of nonadiabatic w a v e pac k et in terferometry using classical tra jectories. Journal of Chemic al Physics 112(17):7345{7354, 2000. [78] Stev e Nielsen, Ra ymond Kapral, and Gio v anni Ciccotti. Mixed quan tumclassical surface hopping dynamics. Journal of Chemic al Physics 112(15):6543{6553, 2000. [79] Ch unCheng W an and Jerem y Sc hoeld. Mixed quan tumclassical molecular dynamics: Asp ects of the m ultithreads algorithm. Journal of Chemic al Physics 113(17):7047{7054, 2000.
PAGE 178
164 [80] Ch unCheng W an and Jerem y Sc hoeld. Exact and asymptotic solutions of the mixed quan tumclassical Liouville equation. Journal of Chemic al Physics 112(10):4447{4459, 2000. [81] Ch unCheng W an and Jerem y Sc hoeld. Solutions of mixed quan tumclassical dynamics in m ultiple dimensions using classical tra jectories. Journal of Chemic al physics 116(2):494{506, 2002. [82] Illia Horenk o, Christian Salzmann, Burkhard Sc hmidt, and Christof Sc h utte. Quan tumclassical Liouville approac h to molecular dynamics: Surface hopping Gaussian phasespace pac k ets. Journal of Chemic al Physics 117(24):11075{ 11088, 2002. [83] D. A. Mic ha and B. Thorndyk e. Dissipativ e dynamics in man yatom dynamics: A densit y matrix treatmen t. International Journal of Quantum Chemistry 90(2):759{771, 2001. [84] Jo el M. Cohen and Da vid A. Mic ha. Electronically diabatic atomatom collisions: A selfconsisten t eik onal appro ximation. Journal of Chemic al Physics 97(2):1038{1052, 1992. [85] ShinIc hi Sa w ada and Horia Metiu. A Gaussian w a v e pac k et metho d for studying time dep enden t quan tum mec hanics in a curv e crossing system: Lo w energy motion, tunneling, and thermal dissipation. Journal of Chemic al Physics 84(11):6293{6311, June 1986. [86] V olk er Engel and Horia Metiu. A quan tum mec hanical study of predisso ciation dynamics of NaI excited b y a fem tosecond laser pulse. Journal of Chemic al Physics 90(11):6116{6128, 1989. [87] J. P Visticot, P de Pujo, J. M. Mestdaugh, A. Lallemen t, J. Berlande, O. Sublemon tier, P Meynadier, and J. Cuv ellier. Exp erimen t v ersus molecular dynamics sim ulation: Sp ectroscop y of Ba(Ar) n clusters. Journal of Chemic al Physics 100(1):158{164, 1994. [88] Thomas Sc hr oder, Reinhard Sc hink e, Suy an Liu, Zlatk o Bacic, and Jules W. Mosk o witz. Photo disso ciation of HF in Ar n HF ( n = 1 14 ; 54) v an der Waals clusters: Eects of the solv en t cluster size on the solute fragmen tation dynamics. Journal of Chemic al Physics 103(21):9228{9241, 1995. [89] M. A. Osb orne, M. A. Ga v eau, C. Gee, O. Sublemon tier, and J. M. Mestdagh. Dynamics of the deactiv ation and desorption of Ba atoms from Ar clusters. Journal of Chemic al Physics 106(4):1449{1462, 1997. [90] E. Czuc ha j, F. Reb en trost, H. Stoll, and H. Preuss. Pseudop oten tial calculations for the p oten tial energies of LiHe and BaHe. Chemic al Physics 196:37{46, 1995.
PAGE 179
165 [91] E. Czuc ha j, F. Reb en trost, H. Stoll, and H. Preuss. Semilo cal pseudop otential calculations for the adiabatic p oten tials of alk alineon systems. Chemic al Physics 136:79{94, 1989. [92] J. P ascale. Use of l dep enden t pseudop oten tials in the study of alk alimetalHe systems. the adiabatic molecular p oten tials. Physic al R eview A 28(2):632{643, 1983. [93] T. Grycuk, W. Behmen burg, and V. Staemmler. Quan tum calculation of the excitation sp ectra of Li He probing in teraction p oten tials and dip ole momen ts. Journal of Physics B 34, 2001. [94] W. Behmen burg, A. Mak onnen, A. Kaiser, F. Reb en trost, V. Staemmler, M. Jungen, G. P eac h, A. Devdariani, S. Tserk o vn yi, A. Zagrebin, and E. Cruc ha j. Optical transitions in excited alk ali + raregas collision molecules and related in teratomic p oten tials: Li + He. Journal of Physics B 29:3891{ 3910, 1996. [95] V. S. Filino v, Y u. V. Medv edev, and V. L. Kamskyi. Quan tum dynamics and Wigner represen tation of quan tum mec hanics. Mole cular Physics 85(4):711{ 726, 1995. [96] B. R. McQuarrie, T. A. Osb orn, and G. C. T abisz. Semiclassical Mo y al quantum mec hanics for atomic systems. Physic al R eview A 58(4):2944{2961, 1998. [97] E. Wigner. On the quan tum correction for thermo dynamic equilibrium. Physic al R eview 40, 1932. [98] K. Singer and W. Smith. Quan tum dynamics and the WignerLiouville equation. Chemic al Physics L etters 167(4):298{304, 1990. [99] Sarah John and E. A. Remler. Solution of the quan tum Liouville equation as a sto c hastic pro cess. A nnals of Physics 180:152{165, 1987. [100] Sarah John and John W. Wilson. Quan tum dynamics as a sto c hastic pro cess. Physic al R eview E 49(1):145{156, 1994. [101] Ra ymond Kapral and Gio v anni Ciccotti. Mixed quan tumclassical dynamics. Journal of Chemic al Physics 110(18):8919{8929, 1999. [102] J. M. Thijssen. Computational Physics Cam bridge Univ ersit y Press, Cambridge, United Kingdom, 1999. [103] P aul A. Fish wic k. Simulation Mo del Design and Exe cution Pren tice Hall, Upp er Saddle Riv er, New Jersey 1995. [104] Mic hael Baer and Rob ert Engelman. A study of the diabatic electronic represen tation within the BornOpp enheimer appro ximation. Mole cular Physics 75(2):293{303, 1992.
PAGE 180
166 [105] T. P ac her, L. S. Cederbaum, and H. K opp el. Appro ximately diabatic states from blo c k diagonalization of the electronic Hamiltonian. Journal of Chemic al Physics 89(12):7367{7381, 1988. [106] A. Thiel and H. K opp el. Prop osal and n umerical test of a simple diabatization sc heme. Journal of Chemic al Physics 110(19):9371{9383, 1999. [107] M. D. F eit, J. A. Flec k, and A. Steiger. Solution of the Sc hr odinger equation b y a sp ectral metho d. Journal of Computational Physics 47:412{433, 1982. [108] C. Leforestier, R. H. Bisseling, C. Cerjan, M. D. F eit, R. F riesner, A. Guldb erg, A. Hammeric h, G. Jolicard, W. Karrlein, H.D. Mey er, N. Lipkin, O. Roncero, and R. Koslo. A comparison of dieren t propagation sc hemes for the time dep enden t Sc hr odinger equation. Journal of Computational Physics 94:59{80, 1991. [109] Andr es Rey es. Density Matrix The ory and Computational Asp e cts of A tomic Col lisions Including SpinOrbit R e c oupling PhD thesis, Univ ersit y of Florida, Gainesville, FL, 2003. [110] Ronald A. Aziz, M. J. Slaman, A. Koide, A. R. Allnatt, and William J. Meath. Exc hangeCoulom b p oten tial energy curv es for HeHe, and related ph ysical prop erties. Mole cular Physics 77(2):321{337, 1992. [111] Da vid J. Griths. Intr o duction to Ele ctr o dynamics Pren ticeHall, Englew o o d Clis, New Jersey 2nd edition, 1989. [112] S. Obara and A. Saik a. Ecien t recursiv e computation of molecular in tegrals o v er Cartesian Gaussian functions. Journal of Chemic al Physics 84(7):3963{ 3974, 1986. [113] S. Obara and A. Saik a. General recurrence form ulas for molecular in tegrals o v er Cartesian Gaussian functions. Journal of Chemic al Physics 89(3):1540{ 1559, 1988. [114] L. E. McMurc hie and E. R. Da vidson. Calculation of in tegrals o v er ab initio pseudop oten tials. Journal of Computational Physics 44:289, 1981. [115] P Sc h w erdtfeger and H. Silb erbac h. Multicen ter in tegrals o v er longrange op erators using cartesian Gaussian functions. Physic al R eview A 37:2834, 1988. [116] J. A. North b y Exp erimen tal studies of helium droplets. Journal of Chemic al Physics 115(22):10065{10077, 2001. [117] Martin H. M user and Erik Luijten. On quan tum eects near the liquidv ap or transition in helium. Journal of Chemic al Physics 116(4):1621{1628, 2002.
PAGE 181
167 [118] D. M. Cep erley P ath in tegrals in the theory of condensed helium. R eviews of Mo dern Physics 67(2):279{355, 1995. [119] K. Birgitta Whaley Structure and dynamics of quan tum clusters. International R eviews in Physic al Chemistry 13(1):41{84, 1994. [120] R. N. Barnett and K. B. Whaley V ariational and diusion Mon te Carlo tec hniques for quan tum clusters. Physic al R eview A 47(5):4082{4098, 1993. [121] F. Stienk emeier, F. Meier, A. H agele, H. O. Lutz, E. Sc hreib er, C. P Sc h ultz, and I. V. Hertel. Coherence and relaxation in p otassiumdop ed helium droplets studied b y fem tosecond pumpprob e sp ectroscop y Physic al R eview L etters 83(12):2320{2323. [122] James Reho, Carlo Callegari, John Higgins, W olfgang E. Ernst, Kevin K. Lehmann, and Giancin to Scoles. Spinorbit eects in the formation of the nahe excimer on the surface of He clusters. F ar aday Discussions (108):161{ 174, 1997. [123] B. Thorndyk e. Cauldr on r efer enc e manual (unpublishe d) 2004. [124] B. Thorndyk e. Qualdr on r efer enc e manual (unpublishe d) 2004.
PAGE 182
BIOGRAPHICAL SKETCH Brian Thorndyk e w as b orn on Octob er 17 th 1969, in Calgary Alb erta, Canada. As he grew up, his mother, Carol Thorndyk e, w as an elemen tary sc ho ol teac her while his father, Gerry Thorndyk e, taugh t highsc ho ol biology and general science. Their p ositiv e academic inruence w as notable, as Brian ended up graduating high sc ho ol with an adv anced In ternational Baccalaureate diploma, and indeed w as oered a full y ear's credit at the Univ ersit y of T oron to in the Departmen t of Ph ysics. Rather than mo v e to T oron to, Brian c hose to remain an additional y ear in Calgary immerse himself in v arious F renc h programs at the Univ ersit y of Calgary and then en ter in to a bac helor's program in ph ysics at an allF renc h univ ersit y in Mon treal. He remained in Mon treal for b oth his bac helor's and master's degrees. After his time in Mon treal, he sp en t a y ear at the Univ ersit y of British Colum bia, w orking on pro jects in the Computer Science and Electrical Engineering Departmen ts. He w ould ha v e remained in computer science in V ancouv er had he not b een en ticed to mo v e to Florida, and complete a master's degree in computer science with Dr. P aul Fish wic k. Tw o y ears later, ho w ev er, Brian decided he had to return to his rst lo v e, ph ysics, and transferred to the Ph ysics Departmen t in the Quan tum Theory Pro ject with Dr. Da vid Mic ha. Tw o y ears b efore graduation, Brian's father w as diagnosed with terminal cancer, and Brian to ok a y ear of semilea v e to sp end time with his paren ts and supp ort them emotionally When his father passed a w a y in Jan uary of 2003, Brian returned to nish his Ph.D., and has since accepted a p ostdo ctorate p osition in the Departmen t of Radiation Oncology Division of Radiation Ph ysics at Stanford Univ ersit y 168



Full Text  
QUANTUM DYNAMICS OF FINITE ATOMIC AND MOLECULAR SYSTEMS THROUGH DENSITY MATRIX METHODS By BRIAN THORNDYKE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 Copyright 2004 by Brian Thorndyke To my father, Gerry Thorndyke ACKNOWLEDGMENTS First of all, I would like to thank my parents, Gerry and Carol Thorndyke, whose unwavering support and love my entire life have allowed me to pursue my dreams. On the academic side, I would like to thank my advisor, Dr. David Micha, for his excellent guidance and encouragement throughout my doctoral work. I would also like to thank the following colleagues in Quantum Theory Project and the Physics department for their friendship and insights during my stay in Gainesville: Jim Cooney, Herbert DaCosta, Alex Pacheco, Dave Red, Andres Reyes, Akbar Salam, Alberto Santana and Zhigang Yi. On a personal level beyond the Physics department, I would be remiss if I didn't express my love and gratitude to Natasha Lepor6. She has been my ..11i sister" for over a decade, and I hope our lives will continue to run with fascinating parallels and intertwine for many decades to come! I'd also like to recognize Albert Vernon who, since our early days in the Computer Science department, has been my partner in our relentless pursuit of aesthetically pleasing code. His friendship has both contributed to some of the best of times and helped me through some of the worst over the last 8 years. Finally, I'd like to express my appreciation to Mike Kuban and Rob Thorndyke for being with me in spirit throughout my Ph.D., and for always providing wonderful avenues of escape and adventure year after year! TABLE OF CONTENTS Page ACKNOWLEDGMENTS ................... ....... iv LIST OF TABLES ........... ...... ........ ...... ix LIST OF FIGURES ................................ x ABSTRACT ................... ............. xiii CHAPTER 1 INTRODUCTION .................. ........ 1 1.1 Overview of Classical and Quantum Dynamics ........... 2 1.2 Approximations to Quantum Dynamics ...... ........ 5 1.2.1 WavefunctionBased Approaches ....... ..... ... 5 1.2.2 Density Operator Approaches ...... .......... 7 1.3 QuantumClassical Liouville Equation ... . 8 1.4 Our Approach ....... ......... ... ...... 10 1.5 Simple OneDimensional TwoState Models . ... 11 1.6 LithiumHelium Clusters ................ .... .. 11 1.7 Outline of the Dissertation ................. .. 12 2 QUANTUMCLASSICAL LIOUVILLE EQUATION: FORMULATION 14 2.1 Introduction .................. ........... .. 14 2.2 Quantum Liouville Equation ............ .. .. 14 2.3 W igner Representation .................. ... 15 2.4 QuantumClassical Liouville Equation . 18 2.5 Effective Potential .................. ........ .. 20 3 QUANTUM CLASSICAL LIOUVILLE EQUATION: COMPUTATIONAL ASPECTS ................ .. ........ ..... 22 3.1 Introduction .................. ........... .. 22 3.2 Trajectory Solution .................. .. .... .. .. 22 3.3 Electronic Basis Set .................. ..... 23 3.4 Nuclear Phase Space Grid ................... 24 3.5 RelaxandDrive Algorithm ................. .. 25 3.5.1 W Independent of Time .... .. 26 3.5.2 W Dependent on Time ....... . ... 26 3.5.3 Velocity Verlet for the Classical Evolution ... 31 3.5.4 Algorithm Details .. ................... 3.6 Computing Observables .. .................... 3.6.1 Operators in an Orthonormal Basis .. .......... 3.6.2 Population Analysis .. ................... 3.6.3 Expectation Values .. .................. 3.6.4 Hamiltonian Eigenstates and Eigenvalues .. ........ 3.7 Programming Details .. ..................... 3.7.1 Orthogonality of Code Development .. ........... 3.7.2 Extensibility . . . . 4 ONEDIMENSIONAL TWOSTATE MODELS .. ........... 4.1 Introduction . . . . 4.2 EffectivePotential QCLE in the Diabatic Representation ..... 4.3 NearResonant Electron Transfer Between an Alkali Metal Surface ..... .............. 4.3.1 M odel Details ... .. .. .. .. .. 4.3.2 Properties of Interest ......... 4.3.3 Results. . . 4.4 Binary Collision Involving Two Avoided Crossings . 4.4.1 M odel Details ............. 4.4.2 Properties of Interest ......... 4.4.3 Results... . . 4.5 Photoinduced Dissociation of a Diatomic System . 4.5.1 M odel Details ... .. .. .. .. .. 4.5.2 Properties of Interest ......... 4.5.3 Results...... . . 4.6 Comparison Using Variable and Constant Timesteps . 4.7 C conclusion . . . Atom and . 41 . 41 . 42 . 47 . 59 . 59 . 60 . 62 . 63 . 63 . 68 . 70 . 72 . 72 5 ALKALI ATOMRARE GAS CLUSTERS: GENERAL FORMULATION 80 Introduction . . . Physical System .. ... .. .. .. .. .. Properties of Interest ........... Hamiltonian for AlkaliRare Gas Pairs ... Hamiltonian for the AlkaliRare Gas Cluster . Electronic Spectral Calculations ...... Electronic Basis of Gaussian Atomic Functions . 5.7.1 Equations of Motion ......... 5.7.2 Overlap Matrix Elements ...... 5.7.3 Kinetic Energy Matrix Elements .. 5.7.4 Coulomb Matrix Elements ...... 5.7.5 Momentum Coupling Matrix Elements . 5.7.6 Dipole Matrix Elements ....... 5.7.7 Pseudopotential Matrix Elements .. . 80 . 80 . 81 . 81 . 84 . 86 . 87 . 87 . 89 . 90 . 90 . 91 . 91 . 92 5.8 Computing the Quasiclassical Trajectory . 5.9 Computational Details ........... 5.10 C conclusion . . . 6 LITHIUMHELIUM CLUSTERS ......... 6.1 Introduction . . . 6.2 Description of the System .......... 6.3 Properties to be Investigated ........ 6.4 Preparation of LithiumHelium Clusters .. 6.4.1 Bulk Helium ... .. .. .. .. .. 6.4.2 Liquid Helium Droplets ....... 6.4.3 LithiumHelium Interactions .... 6.5 Results: Lithium Inside the Helium Cluster . 6.6 Results: Lithium on the Helium Cluster Surface . 6.6.1 Dynamics of Li(2pa) .. 6.6.2 Dynamics of Li(2p7) .. 6.7 C conclusion . . . . 93 . 93 . 97 . 99 . 99 . 99 . 101 . 102 . 102 . 109 . 109 . 119 . 121 . 124 . 130 . 137 7 CONCLUSION .. ............. Effective Potential QuantumClassical Liouville Equation 140 OneDimensional TwoState Models ... . 141 AlkaliRare Gas Clusters .... . 142 Software Development .................. ...... 144 Future Work ....... ...................... 145 APPENDIX A THE CAULDRON PROGRAM .................. ..... 146 A.1 Overview ........................ ....... 146 A.2 Component Descriptions ................ .... 147 A.2.1 Read Input File ... . ..... 147 A.2.2 System: Get Differential Equation Coefficients 147 A.2.3 Propagation: Evolve Single Timestep . ... 149 A.2.4 Properties: Output Properties . 149 A.3 Subroutine Details ................... 150 B SPLIT OPERATORFAST FOURIER TRANSFORM METHOD . C THE QUALDRON PROGRAM. . ..... C 1 O verview . . .. .. C.2 Component Descriptions . ..... C.2.1 Read Input File . . C.2.2 System: Get Hamiltonian Matrix Elements . C.2.3 Propagation: Evolve Single Timestep . . . 151 . 154 . 154 . 155 . 155 . 155 . 155 C.2.4 Properties: Output Properties ... . 157 C.3 Subroutine Details .................. ..... 157 REFERENCE LIST .................. ............. 158 BIOGRAPHICAL SKETCH .................. ......... 168 viii LIST OF TABLES Table page 41 Parameters used in the N.. n f. .:e and Lisurface models. ...... ..42 42 Parameters used in the dual avoided crossing collision model. 60 43 Parameters used in the Nal complex model. ............ ..68 51 Pseudopotential rotation for dfunction mixing. ..... 98 61 Parameters for the HeHe interaction from Aziz (VA). ... 106 62 Parameters for the correction to the HeHe interaction (V) 106 63 Parameters for the eLi interaction. ................. ..113 64 Parameters for the eHe interaction. ................ ..113 65 Parameters for the LiHe core interaction. .. . ..... 114 LIST OF FIGURES Figure page 41 Potential curves for Hamiltonian I: Na incident upon a metal surface. 43 42 Potential curves for Hamiltonian II: Li incident upon a metal surface. 44 43 T(R) at t = 0 au, for the N..ii if.,.e model. This wavefunction is evolved through the SOFFT algorithm. ............. ..48 44 F11 at t = 0 au, for the N..II. f.i.:e model. This PWTDM is evolved through the EPQCLE method. ................. 49 45 T(R) at t = 14000 au, for the N..11. f.,.!e model. . 50 46 Phase space grid points at t = 14000 au, for the N..iii f.,.e model. 51 47 Na populations qi and 7q2 vs. time. ................. 52 48 Li populations mq and 7q2 vs. time. ................. 53 49 Coherence described by Re(iq12) vs. time, for the .. iii f.i:e system. 54 410 Coherence described by Re(M12) vs. time, for the Lisurface system. 55 411 Expectation of position and dispersion for the N.. ,i f.,.:e system. .56 412 Expectation of momentum and dispersion for the N..i if.i.:e system. 57 413 Density function p(R) for the N.. ii f.,:e system. . ... 58 414 Potential curves for the dual avoided crossing collision. ... 61 415 Populations piq and 7q2 vs. time for the dual crossing collision model. 64 416 Coherence described by Re(r12) vs. time, for the dual crossing collision m odel. ........... ...... ........ ..... 65 417 Grid deformation at t = 1400 au, for the dual crossing collision model. 66 418 Probability of transmission in the ground state, for the dual crossing collision model. .................. ..... 67 419 Potential curves for the Nal complex. ................. 69 420 Ionic and neutral populations over time, for the Nal complex. 73 421 Expectation of position and its deviance, for the Nal complex. 74 422 Coherence as a function of time, for the Nal complex. ... 75 423 Phase space grid at the end of the simulation, for the Nal complex. 76 424 Number of steps required by the relaxanddrive algorithm, compared to an estimated number required for a fixed timestep version. 77 61 Schematic of Li(2p) above a He surface. A) Li(2pw). B) Li(2p(r). 101 62 Radial distribution functions for bulk liquid helium .... 106 63 Comparison of the Aziz potential with the effective form. ...... .107 64 Effective HeHe potential. ........ . ..... 108 65 Constraining potential used to keep He atoms from evaporating. 110 66 Temperature fluctuations of the He droplet over time .... 111 67 Helium density profile from the centerofmass of the cluster. 112 68 Adiabatic energy for Li and He as a function of internuclear distance. 114 69 Adiabatic energies for Li and one or more He along the zaxis. 116 610 Adiabatic energy for Li and one or more He along the yaxis. 117 611 Adiabatic energy for Li and a surface of He atoms parallel to the xy plane. . . .. .. .......118 612 Evolution of ground state Li embedded in the center of a He cluster. A) Initial time t = 0 au. B) Final time t = 10,000 au. ...... .120 613 Comparison of Li and He motion within a He cluster. The time scale has been reduced by a factor of 100 for the He curve. ..... ..122 614 Electronic population of Li as it emerges from the He cluster. 123 615 Evolution of Li(2pa) as it recedes from the He cluster surface. A) Initial time t = 0 au. B) Final time t = 33, 000 au. ... 125 616 Mixing of the Li(2pa) and Li(2pr) states at distances where Li(2p) is triply degenerate. ................ ........ 126 617 Electronic population of Li with and without a perturbing electro magnetic field, resonant to the D line. . 128 618 Dipole emission spectra of Li(2pa) as it recedes from the He cluster surface. .... .. .. .... ..... ...... 129 619 Snapshot of Li(2pr) as it interacts with the He cluster surface. A) Initial time t = 0 au. B) Final time t = 67, 000 au. ... 131 620 Electronic population of Li(2p7) as it interacts with the He cluster surface. . .. .. .. .. 132 621 Dipole emission spectrum of Li(2pr) during the first 3000 au. .... 134 622 Dipole emission spectrum of Li(2pr) during the final 3000 au. .... 135 623 Adiabatic curves of Li surrounded by a cubic lattice of He atoms. The parameter R refers to the halflength of the lattice edge. 136 624 Decay of Li(2pr) surrounded by surface He atoms, induced by an EM field with frequency resonant to the Li(2pr < 2sra) transition. .138 A1 Flowchart describing the cauldron program. . 148 C1 Flowchart describing the qualdron program. . 156 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy QUANTUM DYNAMICS OF FINITE ATOMIC AND MOLECULAR SYSTEMS THROUGH DENSITY MATRIX METHODS By Brian Thorndyke M.I. 2004 Chair: David A. Micha M., i r Department: Physics We develop a mixed quantumclassical formulation to describe the dynamics of few and hm..iii body atomic systems by applying a partial Wigner transform over the quantum Liouville equation of motion. In this approach, the density operator becomes a function in quasiclassical phase space, while remaining an operator over a subset of quantal variables. By taking appropriate limits and introducing an effective potential, we derive equations of motion describing quasiclassical nuclear trajectories coupled to quantal electronic evolution. We also introduce a variable timestep procedure to account for the disparity between slow nuclear motion and fast electronic fluctuations. Our mixed quantumclassical method is applied to the study of three simple onedimensional twostate models. The first model represents the photoinduced des orption of an alkali atom from a metal surface, where nearresonant electron transfer is important. A second model explores a binary collision under conditions where two avoided crossings are present. The third model follows the photoinduced dissocia tion of the sodium iodide complex, whose longrange attractive surface results in oscillations of internuclear distance. Quantities such as state populations and quan tum coherence are computed, and found to be in excellent agreement with precise quantal results obtained through fast Fourier transform grid methods. Having validated our approach, we turn to the study of alkali atoms embedded in rare gas clusters, treating the alkali atomrare gas interactions with idependent semilocal pseudopotentials. Light emission from the electronic motion of the alkali atom is derived in the semiclassical limit, and computational methods to render the simulation feasible for a manyatom cluster are discussed. The formalism is applied to lithium atoms in helium clusters, where the cluster configuration and the electronic population dynamics of the lithium atom are monitored over time. We study both the ground and first excited states of lithium, and introduce a resonant electromagnetic field to induce electronic transitions. Our results correlate well with other experimental and theoretical studies on doped helium droplets, and provide insight into the dynamics of an excited lithium atom near a helium cluster surface. CHAPTER 1 INTRODUCTION This work is part of a broader effort to bring new insights into the time depen dence of few and manybody molecular systems. In order to study these systems, we are particularly interested in mixed quantumclassical methods. There are many approaches to combining classical and quantum mechanics,1'2'3 but the underlying theme is to use classical mechanics where quantal descriptions are not essential to the dynamics of the system. By doing so, we can save a tremendous amount of computational time with hopefully minimal expense in accuracy. In our study, we focus on the density operator treatment, which allows for a general introduction of semiclassical and classical limits for some degrees of freedom.4 In the applications being considered, electronic variables are described quantally and the nuclear variables are propagated quasiclassically, but our theoretical treatment of quantumclassical coupling is more general. We are able to validate our methods through extensive study of small test models. These models can be evolved through a fully quantal propagation using fast Fourier grid methods, permitting a rigorous assessment of the accuracy of our approach. As a realistic application, we turn to the study of ground and excited alkali atoms embedded in clusters of rare gas atoms.5'6'7'8'9 Rare gas clusters provide an interesting bridge between fewatom systems and bulk matter. A mixed quantum classical approach allows us to follow the nuclear motion and population state dynamics as the alkali atom interacts with the cluster. By tracking the electronic and induced dipole, we are able to compute the electronic spectra of the alkali atom; and by explicitly introducing an electromagnetic field, we are able to induce electronic transitions. While our formalism is applicable to arbitrary alkalirare gas combinations, we have concentrated on the dynamics of lithium embedded in helium clusters, for which there has been a surge of recent experimental and theoretical activity. 10,11,12,13,9,14 Stable helium clusters doped with lithium atoms have been produced at ultralow temperatures, where the ground lithium atom has been shown to preferentially reside on the surface of the helium cluster. Furthermore, the behavior of the lithium atom subsequent to electronic excitation depends heavily on the orientation of the excited state. Our mixed quantumclassical approach corroborates these findings, and leads to additional insight into the dynamics of these interactions. 1.1 Overview of Classical and Quantum Dynamics The vast majority of chemical and biological processes can be described, in principle, by nonrelativistic quantum mechanics. Within this context, the state of a system of nuclei and electrons is represented entirely by a wavefunction or density operator.15 A Hamiltonian operator describes all interactions between the particles, and can be extended to include environmental components (for example, a boundary or electromagnetic field). When the Hamiltonian does not depend on time, its eigenstates are stationary (up to a phase), and represent timeinvariant configurations of the system. When the Hamiltonian contains time dependence, or the initial state is nonstationary, the molecular system evolves through the action of the Hamiltonian. The solutions to the time independent Schr6dinger equation (TISE) are the eigenstates of the Hamiltonian at any given time. The corresponding eigenvalues are the energy levels available to the system. While conceptually compact, the ana lytical solution of the TISE is not possible for more than two particles in the general case. An enormous body of computational work in chemical and molecular ]1li,i, , is devoted to the numerical solution of the TISE,16,17,18 and accurate ground state energies have been computed for molecules involving hundreds of atoms.19 Knowl edge of the full spectrum of eigenstates would allow one to follow the dynamics of the system as well, but unfortunately it is very difficult to obtain accurate results for excited states, and in any case, the number of states required for accurate computa tions would be prohibitively expensive for most dynamics problems. An alternative is to follow the dynamics directly. Quantum dynamics (QD) follows the evolution of a system in real time. The wavefunction evolves according to the time dependent Schridinger equation (TDSE),20 while the density operator (DOp) evolves through the quantum Liouville equation (QLE).21 These formalisms are equivalent, although statistical ensembles are more naturally described by the density operator. Quantum dynamical calculations are important in understanding the pathways between initial and final states, for exam ple in chemical reactions or molecular collisions, and have become an important complement to modern experiments that use pico or femtosecond light pulses to probe ultrafast dynamics.22 The computational complexity of numerical solutions to full QD, however, severely limits the number of degrees of freedom that can be studied, and incorporation of classical concepts is necessary for most realistic applications. Full QD solutions can be readily implemented by discretizing the wavefunc tion along a multidimensional grid. Since the Hamiltonian contains nonlocal opera tors like kinetic energy, methods such as finite dill. n 1,25,26 or Fourier trans form27,28 are needed to evaluate the action of the Hamiltonian on the wavefunction. In the general case, these techniques require a dense grid to obtain accurate results. Furthermore, their nonlocal nature can result in significant shifts of the wavefunction over time, so that even if the initial state is spatially compact, a large grid is required to accommodate translation. Sparse grid methods29 can alleviate some of the com putational burden of using large grids, and dynamically changing grids30 can better follow the distortion and translation of the wavefunction over time. Alternatively, discretization on a basis set instead of a grid can simplify derivative calculations.31 Ultimately, however, the exponential scaling of the number of grid points or basis functions with the system size renders full QD solutions intractable for more than a few degrees of freedom. Classical molecular dynamics (\!1)) can be derived from the QLE in the classical limit (h  0).32 In the context of molecular simulations, the most basic MD treats nuclei as point particles in phase space, and follows trajectories according to the Hamilton equations of motion.33 Internuclear potentials are derived from ab initio, experimental and empirical results, and the forces on the nuclei are obtained by summing over the partial derivatives of the pair potentials. These classical force calculations are now the bottleneck, and straightforward evolution of N classical degrees of freedom has only O(N2) time complexity.34 Tree methods based on octree spatial partitions35 or multiple potential .::p.ini. .I '. 37 reduce this complexity to O(N log N), with a proportionality constant dependent on the desired accuracy of the simulation. Because of its intuitive nature and computational efficiency, MD is routinely used to study molecular systems with 103 106 atoms.38 The problem with MD is its inability to adequately describe quantum mechanical phenomena, such as charge transfer, electronic excitation, tunneling and zeropoint motion. These effects are ubiquitous in chemical and biological processes at thermal or lower energies, and cannot be completely neglected in most cases.3 An attractive alternative is to construct an approximate QD model that ideally combines the accuracy of QD with the computational efficiency and intuitive sim plicity of MD. In general, this can be done by augmenting MD methods to include quantum features, or by simplifying QD models to incorporate classical or quasiclas sical trajectories.1 There are many variations of this theme, and question remains as to which approach is the most suitable under arbitrary conditions. In Section 1.2, we survey some of the most common approximation schemes, differentiating between methods based on the wavefunction and those centered on the density operator. 1.2 Approximations to Quantum Dynamics 1.2.1 WavefunctionBased Approaches Wavefunctionbased approaches are most useful for the propagation of pure states in closed environments, and have been studied theoretically and numeri cally since the dawn of quantum mechanics. We limit our survey to some prin cipal methods that incorporate some classical concepts to solve the TSDE, includ ing time dependent selfconsistent field (TDSCF) calculations,39,40,41'42 Gaussian wavepacket (GWP) propagation,43'44 surface hopping methods,45 electron nuclear dynamics (END),46'47 and path integral methods.48'49 The time dependent selfconsistent field approximation begins with the wave function written as a product of nuclear and electronic wavefunctions, and uses time dependent variational principles to evolve each wavefunction along a potential averaged over all other wavefunctions. When the nuclear variables are expressed in the classical limit, the TDSCF approximation reduces to a set of electronic wave functions evolving over nuclear trajectories, the trajectories propagating accord ing to effective forces from the quantal system. Many mixed quantumclassical schemes retain this meanfield flavor, where classical trajectories propagate along precomputed or simultaneously evolving electronic states. On the other hand, quantum dynamics has a very different character than classical propagation, and mixed quantumclassical methods vary tremendously in their treatment of quantum classical interactions, initial conditions, and measurement of system observables. Gaussian wavepacket propagation, in its original form, takes the system to be a Gaussian wavepacket in nuclear space, that propagates along a single electronic surface. By locally expanding the potential up to harmonic contributions, equations of motion for the Gaussian wavepacket parameters are derived, which lead to shifts and distortions of the Gaussian over time. The Gaussian center follows classical equations of motion, while its distorting shape results from quantal corrections to the classical motion. The GWP method has been extended to describe nonadiabatic processes with the multiple spawning approach,50 where several wavepackets prop agate on multiple electronic surfaces, and proliferate into additional wavepackets at crossing points. Surface hopping begins with nuclear trajectories evolving on one or more elec tronic surfaces, and reproduces nonadiabatic events through statistically based jumps between the surfaces. In its earliest version, these transitions take place near avoided crossings, but later generalizations allow hops to occur at any time during the sim ulation.51 One drawback to the surface hopping approach is the need to rescale velocities after stochastic transitions, in order to conserve energy. Although accu rate state transitions are often achieved, the dynamics are clearly not representative of the true system evolution, except in a statistical sense. Electron nuclear dynamics is an entirely different approach, that derives equa tions of motion by minimizing the MaxwellSchrodinger Lagrangian density and treating the nuclear degrees of freedom as coherent states52 in the semiclassical limit. By expanding the electronic wavefunction in a basis of Slater determinants, nuclear and electronic motion are combined in a computationally accessible scheme. Path integral methods solve the TDSE for the quantum time evolution opera tor, exp(iHt/h), in the coordinate representation.53 Path integral equations are equivalent to full QD, but unfortunately they are computationally intractable in the general case. Harmonic approximations to the intermolecular potentials sub stantially reduce the complexity of the calculations, making it possible to reproduce the dynamics of a quantum system embedded within an harmonic bath.54'55 For example, variations such as the initial value representation,56'57'58 have successfully described the spinboson model of coupled electronic states in a condensed phase environment. While not applicable to arbitrary systems, these methods are quite useful when the system and its interactions can be cast into the appropriate forms. 1.2.2 Density Operator Approaches There are a number of advantages to using the density operator. First, it provides a convenient representation of mixed rather than pure states.59 Second, many systems can be naturally partitioned into a primary i, 1ii IiI of interest and its surrounding bath. By taking the trace over the bath degrees of freedom, the QLE provides a reduced density description of the primary system interacting with a bath.21'59 Finally, it is possible to include terms directly in the QLE that represent energy dissipation into the environment.60 For these reasons, the QLE is often used where initial conditions are specified as statistical averages, or when departure from a full quantum treatment is necessary due to the size of the system. Among mixed quantumclassical solutions to the QLE, an early approach splits the density operator into a product of Gaussian wavepackets, and derives equations of motions for the Gaussian parameters using selfconsistent field approximations.61 In the semiclassical limit, this method has similarities to generalized forms of GWP propagation, but because it is based on the density matrix, it can naturally treat both open and closed systems. Another method is similar in spirit to END, but uses density matrices rather than Slater determinants to represent the electronic state.62 The quantum degrees of freedom are expanded in some basis set, while the classical degrees of freedom are written as coherent states. The Lagrangian is minimized through the Dirac Frenkel variational principle, and one arrives at equations of motion for a reduced density matrix, coupled to semiclassical equations of motion for the coherent state parameters. The density matrix evolution (DME) method63,64'65 divides the system into a quantum and classical space. The quantum 1i,I 1. iIl is expanded over a nonlocal, orthogonal basis set, while the classical coordinates evolve along trajectories accord ing to the HellmannFeynman force. This method has worked well for simple models involving analytically accessible matrix elements, but is not designed for arbitrary I' II with nonlocal and possibly nonorthogonal basis sets. Eikonal methods66'67'68 exploit the limit of small nuclear wavelengths to sepa rate quantal and classical motion. Through this formalism, the electronic density matrix propagates simultaneously with nuclear trajectories, the latter guided by effective forces from the quantum density. Combined with variable timestep meth ods and travelling atomic function basis sets, the eikonal approach is particularly well suited for the study of binary collision problems. A recent branch of QLE methods uses the Wigner transform (WT) to cast the operatorbased QLE into an equation of motion over phase space variables. The classical Liouville equation for the density function emerges in the classical limit, while judicious application of the WT to a subset of the system variables leads to well defined quantumclassical schemes. In Section 1.3, we outline some of the major approaches found in the literature. 1.3 QuantumClassical Liouville Equation It is fascinating that a phase space representation of quantum mechanics can be developed, that is fully equivalent to the Hilbert space representation. By applying the Wigner transform to both sides of the QLE, an equation of motion is derived for a Wigner function, which is a function of classical position and momentum variables. The Wigner function does not evolve through simple classical equations of motion, however, but involves nonlocal operators in phase space.4 One class of solutions exploits the similarity of these equations of motion with hydrodynamics found in classical li, i. Socalled hydrodynamic quantum dynamics was formulated by Bohm in 1952, and recent developments in the field have revitalized interest in phase space trajectory methods for solving QD.69'70'71'72'73 Unfortunately, these methods, like those involving basis set or grid solutions to the TDSE or QLE, suffer from exponential growth with system size. A mixed quantumclassical scheme arises by dividing the quantum system into quantal and quasiclassical variables, the quasiclassical variables associated with ener gies or masses much greater than the quantal variables. Taking the partial Wigner transform (PWT) of the QLE over only the quasiclassical variables produces an equation of motion for a partially transformed Wigner operator (PTWDOp), which is a function of phase space in the quasiclassical variables, but remains a quantal operator in the space of the quantal variables. After some appropriate approxi mations relating to the different masses of the variables, we arrive at a differential operator equation referred to as the quantumclassical Liouville equation (QCLE).74 While the quantal evolution must still be tackled in Hilbert space, the propagation of the quasiclassical phase space acquires a classical character. If the majority of the system can be described with quasiclasical variables, the computation savings may be considerable. One numerical approach represents the PWTDOp in phase space with a fixed number of delta functions.75,76'77 Specifically, the PWTDop is projected on a basis, resulting in a set of diagonal and offdiagonal functions in phase space. Each function is approximated by a set of delta peaks, that evolve along surfaces constructed from the densities and coherences. Nonadiabatic coupling is implemented through a modified surface hopping procedure, where hops occur between trajectories at curve crossings. An alternative solution uses the same delta peak representation, but rather than use stochastic jumps between trajectories, new trajectories are generated at curve crossings.78 A variation on this approach, the multithread method, spawns new trajectory points at every timestep, but combines the newly generated pools of trajectory points into smaller numbers through energy and other conservation considerations.79,80,81 Recent efforts have been made to combine Gaussian wavepackets in phase space (GPPs) with the trajectory solutions along diagonal and offdiagonal surfaces.82 Rather than represent quasiclassical functions through delta peaks, the functions are expanded as a linear combination of GPPs, thereby more effectively covering quasiclassical phase space. Stochastic jumps between surfaces at crossing points reproduce nonadiabatic behavior, and avoid the need to spawn new GPPs. 1.4 Our Approach In our approach, we introduce an effective potential to the QCLE, to provide a simple numerical procedure for solving the equation.83 Our effective potential quantumclassical Liouville equation (EPQCLE) permits solution in terms of full quantal evolution of the quantal variables along quasiclassical trajectories guided by the quantal state. We introduce a quantal basis and a quasiclassical grid to render the equations suitable for numerical propagation, and implement a scheme to efficiently account for the different time scales of quasiclassical and quantal motion. Our approach shares several positive features with other QCLE methods of solu tion. The QCLE rigorously combines quantum and classical motion, and provides 1,1.I'., i.il ways to incorporate higherorder approximations and quantumclassical coupling. The formalism is based on the density operator, making it attractive for incorporating thermodynamical features or dissipative environments. The QCLE also lends itself to trajectorybased numerical solutions that reduce to the classi cal Liouville equation in the classical limit, providing an intuitive computational framework. Our EPQCLE solution differs from other QCLE primarily in its generality and its scalability. The effective potential can take a variety of forms, and while we have found the HellmannFeynman force to provide the best results for our model '1i I.. it is possible that other forms would provide even better results under different circumstances. Furthermore, an electronic basis is only introduced once the EPQCLE has been developed in terms of partial Wigner operators, and no assumptions are made on the form of the basis set. This flexibility has allowed a uniform treatment of problems ranging from onedimensional models expanded in a twostate diabetic basis, to fully threedimensional atomic clusters in a basis of Gaussian atomic functions. Finally, while the majority of QCLE solutions rely on interactions between the propagating trajectories, our use of an effective potential results in completely independent trajectories whose only connection follows from initial conditions. Such independent trajectory evolution maps optimally to certain parallel architectures, and substantially reduces the cost of propagation, that would otherwise be required to compute trajectory interactions at each timestep. 1.5 Simple OneDimensional TwoState Models We test our methods on a set of onedimensional twostate models, which are simple enough to be evaluated precisely through fast Fourier grid methods. Among these methods, we study a model of an alkali atom approaching a metal surface, where nearresonant electronic transfer is important.84s85 Secondly, we consider a system representing a molecular collision with two avoided crossings, where impor tant interference effects arise as one varies the collision energy.79 Finally, we study a system representing the predissociation of the sodium iodide complex, where the longrange attraction of the excited state results in oscillatory nuclear motion.86 1.6 LithiumHelium Clusters Clusters are .: regates of atoms containing between three and a few thousand atoms. The smallest clusters, or microclusters, have between 3 and 10 atoms, and are just small enough that molecular concepts still apply to some degree. The next range of clusters, or small clusters, have between 10 and 100 atoms, a range where molecular concepts break down and the clusters form shells with nonempty interiors. Large clusters, between 100 and 1,000 atoms, provide the final transition from isolated molecules to the bulk condensed state. Rare gas clusters can be probed by doping them with a chromophore and follow ing the chromophore through various laser detection methods.9'87'88'89 Our study used a lithium atom as the dopant in a pure helium cluster. Semilocal idependent pseudopotentials are used to describe the lithiumhelium interactions, because they are found to accurately reproduce adiabatic energy curves for the lithiumhelium pair for s, p and d symmetries.90'91'92'93'94 By introducing the pseudopotential for malism into the EPQCLE, we are able to follow the nuclear configuration over time. In addition, we are able to probe the evolving electronic energy surface by monitor ing the electronic and induced dipole, and computing the resulting dipole emission spectrum. By applying a resonant electromagnetic field, we are able to stimulate the emission of light by excited lithium atoms near the cluster surface. Both the nuclear configuration of the cluster and the spectral properties of its dopants have been under intense investigation over the recent years, and our contribution gives insight into the dynamics of these interactions. 1.7 Outline of the Dissertation Chapter 2 presents the formalism we have used to explore the dynamics of mixed quantumclassical systems. The ideas are presented for the general case of quantal and quasiclassical variables, and it is shown how the EPQCLE naturally emerges from the partial Wigner transform over the quasiclassical variables. Chapter 3 explores the computational aspects of the EPQCLE. In particular, the application of an electronic basis and a nuclear grid render the equations of motion suitable to numerical solution. We also show how observable quantities can be calculated within this framework. Chapter 4 applies our formalism to three different onedimensional twostate ,1 in,. where the results can be compared to exact quantal simulations based on the Schridinger equation and numerically precise grid methods. Along the way, valuable insights into the accuracy and limitations of our methods are obtained. Chapter 5 presents the formalism for general alkalirare gas clusters. We describe the Hamiltonian for general alkalirare gas interactions, and derive matrix elements for a basis set of Gaussian atomic functions. We further discuss the dipole of the cluster and its matrix elements in detail. Finally, we present computational meth ods used to render the numerical propagation of the equations of motion feasible for small to large clusters. Chapter 6 presents specific results for lithiumhelium clusters. We consider the thermodynamic equilibration of the helium cluster in detail, followed by the evolution of the lithium atom through the cluster and near its surface once the cluster has reached equilibrium. We follow the nuclear motion and discuss the interaction of the excited lithium atom near surface helium atoms, in particular with regard to the nuclear configuration and the dipole spectrum. Finally, we introduce a resonant electromagnetic field to stimulate photon emission after the excited lithium atom embeds itself within the surface. Chapter 7 summarizes the main conclusions obtained in this dissertation. Appendix A discusses the program cauldron developed over the course of this dissertation, to simulate the EPQCLE for the test model systems and the full lithiumhelium cluster. Appendix B derives the numerical evolution of the Schridinger equation through the split operator fast Fourier transform method. The solution obtained by this approach is exact to within numerical precision, but is computationally expensive and thus prohibitive for more than a few degrees of freedom. Appendix C presents the qualdron code, developed alongside cauldron to com pare the results of the test models using mixed quantumclassical dynamics with the full quantal treatment. CHAPTER 2 QUANTUMCLASSICAL LIOUVILLE EQUATION: FORMULATION 2.1 Introduction Rather than focus on the wavefunction of a molecular system, we construct its density operator and proceed from there. The DOp is more general than the wave function, in that it naturally describes a statistical ensemble of quantum states. This is particularly useful when the molecular system interacts with its environ ment, and one does not have a complete knowledge of that environment. The DOp also provides a convenient starting point for deriving mixed quantumclassical meth ods, because the classical limit of the DOp is the classical density function. One particular route to a mixed quantumclassical description of a molecular system is through the the Wigner transform. By applying a partial Wigner transform to the DOp and its equation of motion, one obtains a new representation that is a function in the phase space of the Wigner transformed variables, but remains an operator in the remaining variables. If the PWT is judiciously applied to classicallike (or qua siclassical) variables, then approximations can be made to the equations of motion that provide a classical evolution of the quasiclassical variables coupled to a quantal evolution of the quantal variables. In the case of molecular systems, the quasiclas sical variables are typically the nuclear coordinates, while the remaining variables describe the electronic state. In this chapter, the PWT is described in detail, and its application to molecular systems is outlined. 2.2 Quantum Liouville Equation A system of atoms or molecules is represented by nonrelativistic quantum mechanics as a state vector IT) in Hilbert space. This state vector evolves in time according to Schridinger's equation, ih Ot = HI), (2.1) where H is the Hamiltonian of the full system. If, however, we are studying an ensemble of states {li)} with statistical weights {', }, it is convenient to construct the density operator, = IQ') . (2.2) Taking its time derivative and using Eq. 2.1, we find the density operator to evolve in time according to the quantum Liouville equation of motion, ih = [H, F]. (2.3) There are a number of advantages to using the density operator, but in our case the primary advantage is that it leads directly to the classical density function in the limit h  0. We first discuss the full Wigner transform and its application to the QLE, and then show how the partial Wigner transform can be used to derive a mixed quantumclassical representation in the appropriate limits. 2.3 Wigner Representation The Wigner transform provides a phase space representation of the DOp (the Wigner function) and other quantum mechanical operators. It is defined as the Fourier transform of an operator projected on coordinate space,4'95'96'97 Fw(r, p) (2Ih) dzexp(ip z/)(r z/2r + z/2), (2.4) Aw(r,p) = d'zexp(ip z/h)(r /2AIr + /2), (2.5) where A is arbitrary. The integration generates functions that are local in both coordinate and momentum space, which is important for the emergence of classical features in the development of our mixed quantumclassical method. The prefactors are defined differently for the density operator than for other operators, in order to provide convenient parallels between the Wigner function and the classical Liouville density. Classically, the probability distribution function is well defined. For an Nbody v1, in with coordinates r = (ri, r,..., r3N) and moment p = (pl, I_,... ,P3N), the classical Liouville density p = p(r, p) generates the expectation value of any function A A(r, p),32 (A) drdpp(r,p)A(r,p). (2.6) Quantum mechanically, even the notion of phase space is problematic, as Heisen berg's uncert.,iiii' principle prohibits the simultaneous measurement of position and momentum for a given degree of freedom. However, quasidistribution functions like the Wigner function provide an analogous form of the quantum mechanical expectation value. The expectation value of an arbitrary operator is well known,59 (A) Tr(FA), (2.7) which can readily be seen by expanding the trace over eigenstates of A. However, by Wigner transforming both P and A, we arrive at a form for the expectation value similar to the classical case, (A) drdpfw(r,p)Aw(r,p). (2.8) This can be seen by expanding Fw and Aw in Eq. 2.8, I drdpFw(r, p)Aw(r, p) S(2h) Ndrd I dzexp(ip z/h)(r +z/2,r /2) x f dz' exp(ip z'/h)A(r + z'/2, r z'/2) (2h)3N / drdpdzdz' exp[ip (z +z')/h] xF(r + z/2, r z/2)A(r + z'/2, r z'/2) / drdzF(r + z/2, r z/2)A(r z/2, r + z/2). (2.9) Transforming variables q value, I drdpFw( r+z/2, q' = rz/2, we retrieve the quantum expectation r, p)Aw(r,p) = dqdq'F(q, q')A(q', q) = dqdq'(qPfq')(q'A q'q) f dq(qfPA q) Tr(FA). (2.10) Although Fw is not a probability distribution, various properties justify its classification as a quasiprobability function: I(r)2 dpFw(r,p). (p)2 'drFw(r,p). For any function f(r), (f(r)) drdpFw(r,p)f(r). For any function g(p), (g(p)) = fdrdpw(r, p)g(p). One can also derive the WT of products of operators in terms of the WT of the operators themselves.4 For a general operator P = AB, we find that Aw exp(hA /2i)Bw, (2.11) where A is the bidirectional operator, A  (2.12) Op Or Or Op Expanding the commutator directly, we have [A, B]w = (AB)w (BA)w S Awexp(h A/2i)Bw Bw exp( A/2i)Aw. (2.13) 2.4 QuantumClassical Liouville Equation We have seen that the density operator evolves according to the QLE (Eq. 2.3). Rather than study the solution to the QLE, we will examine the evolution of the Wigner function. When the density operator is transformed over all its variables, we arrive at equations of motion that, in the limit that h 0, reduce to the classical Liouville equation. Here we are interested in performing a partial Wigner trans form of the density operator over a subset of variables (those termed quasiclassical), while leaving the remaining variables (those termed quantal) in their original rep resentation. By taking appropriate limits, we derive equations of motion for the quasiclassical variables, coupled to quantal equations of motion for the remainder. To be specific, we divide the degrees of freedom into N quasiclassical variables Q (Qi,Q2, .* QN) and n quantal electronic variables q = (q, q2,..., ). The Wigner transform is performed over quasiclassical variables only, w(R,P) (27h)N / dZexp(iP Z/h)(R Z/21IR+ Z/2), (2.14) Aw(R, P) dZexp(iP Z/h)(R Z/2PR+ Z/2). (2.15) These are distinguished from fully Wigner transformed operators by virtue of remain ing operators in the quantal space. By taking the PWT of both sides of the QLE, we find the equation of motion for the partially Wigner transformed density operator Sw(R, P),98,99, 100 ih Hw exp(hAf/2i)w Fw exp(hA f/2i)Hw, (2.16) where Ac (2.17) OP OR OR OP To proceed further, we approximate Eq. 2.16 to O(h A cl), 1Fw I I S [I^, F] + ({Hw, Fw} {Fw, Hw}), (2.18) Ot ih 2 where {A, B} is the Poisson bracket, {A, B} A A ciB. (2.19) This is an approximation within the quasiclassical space, and reduces quantal motion of the quasiclassical degrees of freedom, but keeps dependence of the quasiclassical variables on the quantal state through the interaction potential. It is an appropriate truncation when the quasiclassical variables are associated with a much greater mass than the quantal variables.75'101'80'79 For a specific example, consider a Hamiltonian for a molecular system composed of kinetic, potential and interaction terms, p2 2 H +V(Q)+ ()+ '(4, + Q), (2.20) 2M 2m where P is the nuclear momentum operator, V(Q) the nuclear potential, p the elec tronic momentum, v(q) the electronic potential, and V'(q, R) the electronicnuclear coupling. By interpreting the nuclear variables Q as quasiclassical and taking the PWT over the nuclear variables, we find the partially Wigner transformed Hamil tonian, 2 + V() + v() +V'(,). (2.21) Defining fH +p + ) + V'(4, R), 2m V V(R) + V'(,R), (2.22) (2.23) we get the quantumclassical Liouville equation of motion for Fw(R, P), OFw Ot [I q p Fw] I+ + (2.24) ih M 2 OOR OP aP aR 2.5 Effective Potential The third term on the RHS of Eq. 2.24 presents a challenging obstacle to solving the equation numerically. One problem common among many proposed schemes is the requirement of computationally demanding algorithms to evaluate the partial derivatives of the PWTDOp. To avoid this problem, we introduce an effective potential V(R, P) in Eq. 2.24, 9Pw 1 P avP w ap 1V OE w [H i, Fw ] + Ot ih M OR OR OP { ( R aR a _R aR (2.25) 2 OR OR OP OP iR ORf v The introduction of V becomes computationally useful insofar as we can neglect the fourth term in Eq. 2.25. While the choice of V is clearly arbitrary, we have found optimal results by setting the expectation value over quantal variables of (aRVaRv) to zero, Trqu 1W ( RV a VI OR OR (2.26) In this way, we retrieve the HellmannFeynman force, 9 Trq, q['wwv/9Rj 0 Trq u[P (2.27) OR Trq. [fW] FHF(R,P). (2.28) The denominator in Eq. 2.27 is itself a function of quasiclassical phase space, which differs from effective potential approaches involving a single normalized density oper ator. Neglecting the fourth term in Eq. 2.25, we find an approximated but compu tationally advantageous effectivepotential QCLE (EPQCLE), 0f1w 1 P 80w 0 aw Ot i[H F M OR F P The approximations used to derive the EPQCLE substantially reduce the quantal character of the quasiclassical solution space. In contrast, the best adia batic methods, based on the BornOppenheimer separation of nuclear and electronic motion, retain full quantum nuclear dynamics along adiabatic curves.102 A princi pal advantage to the EPQCLE is its nonadiabatic character, which is capable of describing nonadiabatic events when the BornOppenheimer limit no longer applies. A further benefit is its suitability for solution using trajectory methods, greatly increasing the size of problems amenable to numerical analysis. The theoretical and computational aspects of the trajectory solution to the EPQCLE are discussed in the next chapter. CHAPTER 3 QUANTUM CLASSICAL LIOUVILLE EQUATION: COMPUTATIONAL ASPECTS 3.1 Introduction The EPQCLE is a partial differential operator equation in the quasiclassical variables and time. One way of solving this kind of problem is to represent the operators as matrices on a large grid, and evolve the matrices using finite difference or spectral methods. The drawback to this approach is that a very dense grid is required for numerical accuracy, and a very large grid is necessary if the quasiclassi cal density shifts location appreciably as it evolves in time. Since the grid dimension varies directly with the classical degrees of freedom, a multiparticle system presents very serious numerical difficulties. Moreover, finite difference and spectral grid solu tions are inherently difficult to parallelize, as substantial communication is required between processors regardless of the division of computational labor. An alterna tive approach, applicable to the class of partial differential equations to which the EPQCLE belongs, is to follow trajectories in phase space as the system evolves. Only the important trajectory points are represented, so that the (moving) grid maintains a minimal size. In this chapter, we explore the trajectory approach and see how the EPQCLE can be solved in an efficient and even parallel manner. 3.2 Trajectory Solution We can formally solve Eq. 2.29 by following trajectories in classical phase space, with R and P becoming functions of time, dR P (3.1) dt M' dP d FHF(R,P). (3.2) dt One could follow any paths in phase space, but by using those II'... . .1 in Eq. 3.2 we are able to transform Eq. 2.29 into an ordinary differential equation in time. Inserting Eq. 3.2 in Eq. 2.29 and moving the partial derivatives to the LHS, we derive the change of the PWTDOp along the quasiclassical trajectories, dF 1 Sfi[HF] (3.3) dtih Note that we have omitted the subscript on the PWTDOp for notational conve nience, and from here onward will continue to label all PWT operators without this subscript. Eq. 3.3 remains a formal solution, and before solving it numerically we must discretize the equations in both quantal and quasiclassical space. This is the subject of the next two sections. 3.3 Electronic Basis Set Let us introduce an arbitrary basis, { tions, a Gaussian basis is used, but for now we consider the basis to be general, and not necessarily orthogonal or normalized. Converting to matrix notation, we let 1I4) be the row matrix, 14) (1 I) 1 ). (3.4) Then we can expand our operators, PF I4) Sl (FIIl4)S1 (4,, (3.5) F A 4) S1 (4 AI4) S I(4 (3.6) A where S is the overlap (44). Projecting Eq. 3.3 on this basis set, and setting h 1, we obtain dT r(iH q _t)S1 Sl(iH' + )F (3.7) dt where we have used the notation, dR dP S (fl ( //&R4) + (4d/&P4). (3.8) dt dt 3.4 Nuclear Phase Space Grid Although we've transformed the EPQCLE into a discrete representation in electronic space and are now dealing with the partially Wigner transformed density matrix (PWTDM) instead of operator, we must still discretize the quasiclassical phase space. To this end, we choose a set of initial grid points {(Ri, Pi)}. Their distribution should approximately cover the domain of F(R, P), and should be suf ficiently dense to well represent the evolution of the PWTDM. In practice, the grid should be adjusted until convergence is achieved. Once a grid is chosen, the grid points follow nuclear trajectories {(Rj(t), Pj(t))} according to Eqs. 3.1 and 3.2: dRj P1 S (3.9) dt M' dP FHF(RJ, P). (3.10) dt At the same time, Eq. 3.7 becomes a set of uncoupled equations, one for each trajectory: dr F(iHq Q)S S 1 (iH + (3.11) dt r^H s H (3.11) where r r((Rj(t),Pj(t)), (3.12) Hj H(Rj(t),Pj(t)), (3.13) Qfj F(R,(t),Pj(t)). (3.14) Each trajectory follows the evolution of the PWTDM along that path, independent of the other trajectories. While one cannot expect coherence between the classical degrees of freedom to be represented by this approach, there are some substantial computational advan tages. In particular, the scheme can be optimally ported to a parallel processor, whereby each processor independently evolves a single trajectory; communication between the processors is unnecessary. 3.5 RelaxandDrive Algorithm Before the trajectory solution can be implemented, a detailed propagation scheme needs to be specified. One would expect that with a sufficiently small timestep At, all propagation methods would converge to the same results, provided roundoff error were not significant. However, the paths to convergence will certainly differ, in that methods with higher accuracy can use larger timesteps. The relax anddrive method, developed originally by Micha and Runge,66'67'68 incorporates the rapid electronic oscillatory behavior with the relatively slowly evolving nuclear variables in an accurate, variable timestep scheme. The relaxanddrive procedure has been shown to give excellent results for a wide variety of twobody collision problems. For completeness, we review its details next. First of all, since all trajectories are propagated analagously, it suffices to look at the evolution of a single trajectory. Accordingly, let us rewrite the EPQCLE for a single trajectory, S =. W(t)F(t) F(t)Wt(t), (3.15) dt where we have defined W S1(H + i+). (3.16) We wish to propagate from initial conditions Wo = W(to), Fo = r(to). If W is independent of time, we find a numerical solution that is exact up to machine precision. If W depends on time, but is slowly varying with respect to the timescale of F, we can propagate to a high degree of accuracy by linearizing Eq. 3.15 in time and incrementing in small timesteps At = tl to. 3.5.1 W Independent of Time When W is independent of time, W(t) = Wo, we can formally solve Eq. 3.15, F(t) Uo(t,to)roUt(t, to), (3.17) where U(t, to) = exp[iWo(tto)]. (3.18) If we diagonalize Wo, Wo = TAT1, (3.19) we can rewrite Eq. 3.17 as, F(t) T[TUo(t, to)T]T o(Tt)1[TtUt(t, to)(Tt)1]Tt. (3.20) The exponential matrices can be formed analytically, since T1Uot, to)T exp[iA(tto)], (3.21) and F(t) can be computed in the time complexity of the diagonalization of Wo. 3.5.2 W Dependent on Time When W depends on time, we separate F into a reference F and correction Q term, r(t) r(t) + Q(t). (3.22) The reference density is propagated by the time independent Wo of Section 3.5.1, dro(t) dt The evolution of Q(t) is formed by inserting Eq. 3.22 into Eq. 3.15, (3.23) .d i(ro + Q) dt (Wo + AW)(r Q) (r + Q)(Wo + AW)t, (3.24) where AW W Wo. (3.25) Using Eq. 3.23 we obtain AWro + WoQ + AWQ FOAWt QWt QAWt. (3.26) Transforming to the local interaction picture, where A UoALU , UoAWLUoUoFrLU +UoAWLUoUoQLUt UoFiUoUoAW Uo UoQLUoUoAWLUo. (3.28) We can simplify Eq. 3.28 by multiplying on the left by U 1 and on the right by (U )1 to obtain DL + [AWLU UoQL QLUoUoAWj], (3.29) where AWLUoUoFL FLUoUoAWL. dQ dt we get iU dQL iUo Uo dt 0 (3.27) dQL dt WoF(t) ro(t)Wt0 DL(t) (3.30) Formally solving for QL, Q(t) (t) + dt'[AWL(t')U(t', to)Uo(t', to)QL(t') to QL(t')Uo(t', to)Uo(t', to)AWt(t', to)], (3.31) where t Qf (t) f dt'[AWL(t')U(t', to)Uo to)F (t') to ro t')U (t', to)Uo(t', to)AW (t', to)]. (3.32) Solving by iteration, Eq. 3.31 becomes t Q t) Qf() f + dt'[AWL(t')Utt',to)Uo(t', to) Qf(t) to Qf (t')U (t', to)Uo(t', to)AW (t', to)] (3.33) Neglecting the second term and higher for small timesteps At = ti to, tl QL(tf) dt'[AWL(t')Ut(t', to)Uo(t1', 10) (tl') to ro (t')Uto(t', to)Uo(t', to) AW (t', to)]. (3.34) Converting back to the original representation, AL = UoA(U)1, (3.35) we have t1 Uo(tl, to)Q(t1)Uo(t1, to)1 dt'[Uo l(t, to)AW(t')ro(t')Ut(t', to)1 to Uo1 (t', to)Fo(t')AW (t')U (t', to)1]. (3.36) Multiplying Eq. 3.36 on the left by Uo(ti, to) and on the right by Ut(ti, to), and noting Uo(t, to)Uo (t', to) exp[iWo(ti exp[iWo(ti Uo(tl, t'), exp[iWV(t' exp[iW(ti  to)] exp[iWo(t' (3.37) U (ti, t'), (3.38) we get an approximate correction term, Q(ti) to dt'Uog(t t')D(t')U (tl, t'), to where D(t') AW(t')F(t') F(t')AWt(t'). We can compute Q(tl) by quadrature if we assume D(t) = D(t1/2), to t < ti where t1/2 to + At/2. We then have Q(ti) tl S dt'Uo(t, t')D(tl/2)U (t t'), to (3.39) (3.40) (3.41) (3.42) Ut (t, to)I u (ti, to) to)] exp [^iWo (ti which we can rewrite as Q(ti) T Idt'exp[iA(ti to t')]DT exp[iAt(ti T D(tl/2)(Tt)1 Examining the elements of Q, t')] T k Sdt'exp[iAo(tl oi f dt' exp[ to A )(ti t')]DL Ttk 1 exp[i(Ai A )At] T (A m A )k where we have used the notation Ai Aii. Reverting back to matrix form, we have Q(ti) TXTt, exp[i(A A )At] 1 T A, (TAm T D(tl/2) )1. The full density matrix at tl is then simply, (3.48) tl dtT[T Uo(t, t')T]DT[TU(Tt)1]T to where t')] Tt, (3.43) Qjk(tl) (3.44) Til Tm Im i Ti1 l ZTn i(Al (3.45) where (3.46) XDm DT (3.47) t')]DI exp[iA (t, F(ti) r(ti) + Q(ti). 3.5.3 Velocity Verlet for the Classical Evolution In order to complete the relaxanddrive algorithm, we need to account explicit for the propagation of the classical variables. We do this by assuming that during each timestep, R and P are advanced by the reference density F, and that the correction term contributes negligibly to the classical propagation. The precise nature of the classical propagation is independent of the relaxanddrive algorithm, although to keep accuracy to O(At2) it is necessary to integrate using an algorithm like velocity Verlet or RungeKutta. We choose to use the velocity Verlet method, which is accurate to O(At2) and is selfstarting. We proceed by advancing the classical positions, P (t) 1 dP(t) 2 R(t + At) = R(t) + At + t2. (3.49) M 2M dt The last term of Eq. 3.49 is the acceleration term, which comes from the effective potential. Having advanced the positions, we calculate the acceleration at the new location, and advance the moment, P(t + At) P(t) + At. (3.50) In the context of the relaxanddrive procedure, the classical coordinates are advanced in two steps of At/2, in order to compute the correction term Q. Variable timestep. As the propagation proceeds, the initial timestep may no longer be appropriate. This often occurs when the system enters a region where the magnitude of potential interactions changes so that the reference density fluctuates at a different rate. An example is the collision of two atoms, where large timesteps can be taken at large distances, but small timesteps are required at close range. We can monitor this fluctuation by observing the correction term Q. As Q becomes too small (large), we need to increase (decrease) the timestep to maintain efficiency (accuracy). To this end, we define a correction measurement, if 2 e i max (3.51) At the end of each timestep, we evaluate e. If c is less than some threshold, say cl, we discard the step and use a new c 3 x c. If c is greater than some threshold, .1,v c, we discard the step and use a new  e/2. By either multiplying by 3 or dividing by 2, we avoid any oscillation (and thus infinite loop) between a pair of timesteps. 3.5.4 Algorithm Details The algorithm can be divided into an outer piece (say, Main), which calls an inner piece (say, Propagate). They are described in point form as follows: Main: To propagate F from tA to tB, given At = Ato and {fl, ej}, 1. At min{At, tB tA}. 2. Propagate, with to = tA. 3. Compute = n.i:: ; Qj/  l2. 4. Is e < I and tA At < tB? YES: Reset variables (matrices and functions return to their values at to), set At  3 x At and return to step 1. 5. Is e > c,? YES: Reset variables to to, set At + At/2 and return to step 1. 6. IstA At YES: Set tA ` tA + At, and return to step 1. Propagate: To propagate F from to to tl = to + At, 1. Initialize variables: Wo = W(to), Fo = ro(to) = r(to), R(to), P(to), At. 2. Diagonalize Wo = TAT. 3. Advance classical variables by half timestep, {R(tl/2), P(tl/2)}, where t1/2 t + At/2, using initial reference density Fo(to) and Wo. Compute W(tl/2) using R(ti/2) and P(tl/2). Compute or(tl/2) Texp[iA(At/2)]Tlro(Tt)1 exp[iAt(At/2)]Tt. Compute Q(ti) by assuming D(t) = D(tl/2), to < t < ti. (a) Calculate AW(tl/2) = W(t/2) Wo. (b) Calculate D(tl/2) AW(tl/2)Fr(tl/2) F(tl/2)AWt(tl/2). (c) Compute DT = T1D(tl/2)(Tt)1 (d) Compute X, where XIm. exp[i(A A )At]  tA DlT A, Amtn (e) Compute Q(ti) = TXTt. 7. Advance classical variables to full timestep, {R(ti), P(ti)}, using Fo(tl/2) and W(ti/2). 8. Compute W(ti) using {R(ti), P(ti)}. 9. Compute Fr(ti) Texp(iAAt)Tlro(Tt)1 exp(iAtAt)Tt. 10. Compute F(ti) F (ti) + Q(ti). 3.6 Computing Observables 3.6.1 Operators in an Orthonormal Basis Up to this point, we have considered the general basis 4) without any condition of orthogonality or normality. We can transform this basis to an orthonormal one, V'), through a L6wdin transformation, I4)S1/2 (3.52) Since quantal traces are naturally formulated in orthogonal bases, it is useful to express the relationship between matrix representations of operators in both bases. For the density operator, r I\'}r'{^' I=)S1/2r'S1/2(4 )  F)r 1. (3.53) Thus we equate S S1/2F'S1/2. (3.54) Since the density matrix has been defined differently than the matrix representations for general operators, we also consider the representations for a general operator A, A ')A'({' = 1)S1/2A'S1/2( I4)S'AS (). (3.55) Thus A S1/2A'S1/2. (3.56) 3.6.2 Population Analysis The population is naturally defined in an orthonormal basis, such that the population of state i is the ith diagonal element of the orthonormal representation of the density matrix, ri [Sii [S1/2rsl/2]ii. (3.57) In the case of the PWT representation we must integrate over nuclear phase space as well, so that Eq. 3.57 becomes S= dRdP[S1/2(R, P)F(R, P)S1/2 R, P) (3.58) 3.6.3 Expectation Values The expectation value of a general operator Aw is found by taking both the quantal and classical trace of the product of the operator with the PWTDM, (Aw) = Tr[AwFw] STrciTrqu[AwFw] dRdPTrqu[Aww]. (3.59) The quantal trace is naturally computed in the orthonormal basis, but as we now show, the nonorthonormal basis representation can also be used: Trqu[AwFw] = Trqu[A' F'] STrqu[S1/2AwS1/2S1/2FrS1/2] STrqu[AwFw]. (3.60) 3.6.4 Hamiltonian Eigenstates and Eigenvalues Using the orthonormal representation of the Hamiltonian, H' = S1/2HS1/2, (3.61) we can compute its eigenvalues by diagonalizing the matrix, L1H'L = H', (3.62) where H' is the diagonalized Hamiltonian and contains the energy eigenvalues along its diagonal. Since the Hamiltonian is Hermitian, its eigenvalues are real, so that LtH'(L')t = H, (3.63) and the diagonalizing matrix L is seen to be unitary, Lt = L. (3.64) The columns of a unitary matrix are orthogonal, and since the columns of L are the eigenstates of the Hamiltonian, we see that the eigenstates produced by diagonalizing the orthonormal representation of the Hamiltonian are also orthonor mal. This is useful in the case of degenerate eigenstates, as they are automatically orthogonal and no additional procedures are needed to ensure orthogonality in the degenerate subspace. 3.7 Programming Details When designing a computational package, it is desirable to ensure the code remains orthogonal and extensible throughout the design and implementation. Mod ern programming languages use objectoriented concepts to achieve these goals,103 but unfortunately a great deal of computational work is built on older procedural languages and cannot be readily incorporated into an objectoriented scheme with out substantial effort. Indeed, as much of this legacy code has been developed over many years and has undergone extensive t. l iii. it is enticing to adhere to the older languages in which they were written and incorporate them directly. In the devel opment of code for this research, a compromise was found by using many advanced features found in Fortran 90, but maintaining a coding style which permit straight forward integration of legacy Fortran 77 code. Details of the package (cauldron) are found in Appendix A, and in the remainder of this chapter we briefly overview the programming principles used in the development of the package. 3.7.1 Orthogonality of Code Development By orthogonality of code development, we mean that different aspects of the code can be developed independently. Thus one may decide to build a completely dif ferent algorithm than relaxanddrive, for example, to propagate the mixed quantum classical system, but be able to do so without changing aspects of the code which define the system, compute its properties, read configuration files, generate output files, and so forth. This helps ensure that once a version of the code works well, changes to its components will be less likely to introduce errors. This aspect of program development is crucial for software designed in a team environment, but also very useful for the solitary designer when the problems and their solutions may rapidly change. Orthogonality has been kept within cauldron through judicious use of variable and subroutine naming conventions, and by building a solid hierarchy of directories and subroutines from the beginning. 3.7.2 Extensibility In scientific work, the systems studied and the solutions used are constantly (1I.ii,il,. as progress is made in understanding the solutions, and new problems arise. One way to help maintain flexibility, which has been used throughout cauldron, is to ensure that systems are represented as generically and as dynamically as pos sible. Generic code attempts to represent the fundamental aspects of all molecular  1. ii, . for example, by the same set of variables and .iili.,. When new com ponents (e.g., new kinds or numbers of nuclei) are added to the system, the same variables are used, and it is only the interpretation of the results that differs from ,i1. I, to system. By dynamic representations, we refer to the defining of variable size at runtime rather than fixing the size at compilation. The major benefit in this comes from being able to implement models of varying sizes without recompiling the code and creating a new executable for each system studied. Systems can then 38 be defined completely within input files, for example, preserving the polished exe cutable without modification. Fortran 90 encourages dynamic memory allocation, which has been used to great advantage in cauldron to provide very extensible code. CHAPTER 4 ONEDIMENSIONAL TWOSTATE MODELS 4.1 Introduction While our ultimate goal is to study realistic threedimensional models of alkali atoms embedded in rare gas clusters, the complexity of these systems places full quantal solutions out of computational reach. On the other hand, simple models can sometimes capture elements of larger and more realistic systems, and provide a rigorous basis for validating approximate numerical methods. In this chapter, we study the dynamics of three simple models involving two electronic states and one nuclear coordinate. The first represents photoinduced desorption of an alkali atom from a metal surface, where nearresonant electron transfer is important. The second models the collision between two nuclei in a framework involving two avoided crossings. The third models the photoinduced dissociation of the Nal complex, where oscillatory motion between neutral and ionic states is observed. Because of the limited size of these models, in each case we are able to propagate a grid solution to the TDSE, and thus compare our EPQCLE approach to the dynamics of the full quantal system. For all three models, we will see that the mixed quantumclassical methods provide qualitatively, and often quantitatively similar results to the full quantal evolution. 4.2 EffectivePotential QCLE in the Diabatic Representation In Chapter 3, we derived the EPQCLE for an arbitrary basis. Here, we consider the specific case where the system is described in an orthonormal diabetic basis. There are many varieties of diabetic bases,104'105 but here we refer to the strictly diabetic representation {()d} where the momentum coupling vanishes,106 ({Qd =/Ol d) = 0. (4.1) We also assume that the basis does not explicitly depend on P or t, so that Q = 0. Since the basis is orthonormal, the overlap is unity, and Eq. 3.7 reduces to the simple form, dF S [Hqu, ]. (4.2) In the diabetic representation, the effective force is also simplified, since the opera tors in the quantal trace can be replaced directly with their matrix representations, Trqu [FOV/OR] FHF (4.3) Trqu [r] For our test models, the partial derivatives of the potential can be calculated ana lytically at each grid point in phase space. Since the PWTDM is propagated along these grid points, the product in the numerator of Eq. 4.3 is computed through matrix multiplication, while the quantal trace is calculated by summing over the diagonal components of this matrix product. The quantal trace in the denominator, on the other hand, is simply the sum over the diagonal components of the PWTDM. Had we used a different basis, we would not have been able to simplify the EPQCLE in this way. However, the real advantage to the diabetic representation is that for very small systems, it lends itself to a fully quantal numerical solution through the propagation of the TDSE on a grid. One scheme, the split operator fast Fourier transform (SOFFT) method, splits the Hamiltonian into its kinetic and potential components, and uses the fast Fourier transform to compute the evo lution due to the kinetic terms.107'108 While this method is very accurate, it is also intractable for systems with more than a few degrees of freedom. However, because our models are simple, the SOFFT procedure provides an excellent test of the accu racy of the EPQCLE. A complete description of the SOFFT is given in Appendix B, and the code implementing the SOFFT (qualdron) is described in Appendix C. 4.3 NearResonant Electron Transfer Between an Alkali Atom and Metal Surface 4.3.1 Model Details In the first of our test systems, we consider a model describing the nearresonant electron transfer between an alkali atom (Ak) and a metal surface (I\) at thermal energies. The model consists of two diabetic surfaces corresponding to a ground state of neutral components Ak + M (state 1) and an excited state for ionic components Ak+ + M (state 2), which cross at short distance. The surfaces and interaction term are given by,85,84 H11(R) = Uo{exp[2a(R Ro)]+ 2exp[a(R Ro)]} /2, (4.4) H2(R) Uo{exp[2a(R Ro)] 2exp[o(R Ro)]} + /2, (4.5) H12(R) = cexp[a2(R 1) 2]. (4.6) Here, R is the distance between the metal surface and the nuclear center of the Ak atom, and is the quasiclassical variable over which we take the PWT. The ionic curve, H22(R), is a Morse potential with a binding energy Uo. The repulsive neutral curve, H11(R), is offset relative to the ionic curve to give an excitation potential Ac. The strength of the coupling term, H12(R), is characterized by 3 and peaks at the crossing R = R, between H11 and H22. The initial state is formed by approximating the ionic surface as a harmonic potential around its minimum R = Ro, and finding the lowest bound vibrational state within that (harmonic) well, (R) = ( ) exp R2 exp[iPo(R Ro)]. (4.7) 7ff2 $U Table 41. Parameters used in the N..U1 f...:e and Lisurface models. Hamiltonian I Hamiltonian II Parameter value (au) value (au) Uo 0.025 0.184 a 0.4 0.4 e 0.005 0.147 Ro 5.0 5.0 R, 12.5 9.0 Po 0.0 0.0 a 0.233153 0.1i f 0.15 0.15 M 42,300 12,800 The PWTDOp is the PWT of Eq. 4.7 over R, giving a Gaussian density in (P,R), F(P,R) I RoR 2( _Po)2 (4.8) At t = 0, the electronic state is promoted by a sudden optical excitation to the repulsive neutral potential, so that the PWTDM becomes l(P, R) = 1 (R Ro 2 U(P PO P (4.9) with F12 F21 = 22 0. The simulation follows the spontaneous decay of this state. The parameters used in the calculation are shown in Table 41, where we con sider two model Hamiltonians: (I) N..II[ f.,:e and (II) Lisurface. The diabetic potentials for Hamiltonians I and II are shown in Figures 41 and 42. 4.3.2 Properties of Interest Populations. We can follow the populations over time by taking the full classical trace over either diagonal element of the PWTDM, r] J = fdRdPri(R, P). (4.10) (7 43 0.02 H22 22 ro  0.015 \H12 0.01 0.005 s 7            H0    " < ^  CT 0 C 0.005 LU 0.01  0.015 0.02  0.025 0 5 10 15 20 25 30 R (au) Figure 41. Potential curves for Hamiltonian I: Na incident upon a metal surface. 44 0.15 H11 \H22 H12 0.1 0.05 5 'i  " ) 0 0.05 0.1 0.15 0 5 10 15 20 R (au) Figure 42. Potential curves for Hamiltonian II: Li incident upon a metal surface. Since the system begins in the repulsive (neutral) state, one expects a certain per centage of the population to fall to the attractive (ionic) state as the system passes through the region of nonnegligible potential interaction. Coherences. The coherence between the neutral and ionic state is described by the real and imaginary components of the offdiagonal terms of the PWTDM, for example i= j dRdPiFj(R, P). (4.11) Coherence is a purely quantum phenomenon, and one measure of the 11,.,lil of a mixed quantumclassical method is the degree to which it maintains coherence. Position expectation values. One observable we can study is the expecta tion of position, (R) =Tr[F(R,P)R] I dRdPTrqu[F(R, P)]R. (4.12) We can also measure the dispersion, aR [((R (R))2)]1/2 [f dRdPTrqu[F(R, P)](R (R))2 (4.13) Momentum expectation values. Similarly, we can compare the expecta tion and dispersion of moment, (P) = Tr[F(R, P)P] f dRdPTrqu[F(R, P)]P, (4.14) p = [((p (P))2)]1/2 / dRdPTrqu[F(R, P)](P (P))2 (4.15) / 1/ Probability density. We can compute the probability density p(R) from the PWTDM by taking the trace over quantum variables and moment, 00 p(R) = dPTrq[F(R, P)]. (4.16) 00 In practice, since the grid in phase space quickly deforms as the system evolves, we must find a way of approximating this integral. One procedure is to determine the support of the PWTDM in quasiclassical phase space, [Rmin, Rmx] x [Pmin, PFia]. We then divide this space into NR x Np equisized bins {bi}, such that bin bij spans the rectangular region, Rn (i 1 Rmax Rmin Rmax Rmin [Rmin + (i 1) Rmin + i ] NR NR x [PTi + (j ) Pmin + J ]. (4.17) Np Np We then assign a value to each bin, pij, which is the weighted sum of all Nij trajec tory points which fall within that bin, Nij r F(Rk, Pk) P 1 (4.18) We can determine the matrix probability density from pai by summing over all bins containing a given position R, p(R) = VpR, VR b. (4.19) Finally, we compute the probability density by taking the quantal trace over the matrix probability density, p(R) = Trqu[p(R)]. (4.20) This probability density can be compared to the density function obtained from the SOFFT simulation, p(R) = I(R) 2 + 2( R)2. (4.21) Wavefunction and PWTDM. Of course, the evolution of the quantum wavefunction can be contrasted directly with the evolution of the PWTDM. How ever, we can also observe the distortion in the phase space grid used by the EPQCLE method. While the grid is initially uniform, it changes shape in an interesting way because of the action of the effective potential. 4.3.3 Results Figures 43 and 44 show the initial wavefunction and its PWT, respectively. Note that the PWTDM formed from a Gaussian wavepacket is a Gaussian function itself, albeit in twodimensional phase space. Figure 45 shows the wavefunction at t = 14000 au, having been propagated through the SOFFT method. The PWTDM evolves in phase space through the EPQCLE, and in Figure 46 we show the grid in phase space at the final time. While substantially distorted from its initial uniform distribution, we notice that the points are globally positioned along a straight line in phase space. This reflects the .ii. ,,,.1, . i' state, where each point is subject to a vanishing HellmannFeynman force, and thus propagates at constant velocity. The observables are presented in Figures 47 to 413. There are a number of interesting things we can glean from these plots. From Figures 47 and 48 we see that as the atom moves away from the metal surface, much of its population shifts from the neutral to ionic state. In the case of Na, approximately 2/3 of the neutral population shifts, while for Li the transfer is total. This may reflect the stronger interaction coupling involved with Li. Figures 49 and 410 describe the coherence between the states. For the N. ii f.i.e system, the coherence remains large for long times, while for the Lisurface 3 3.5 4 4.5 5 5.5 6 6.5 7 R (au) Figure 43. T(R) at t = 0 au, for the N..i i f..e model. This wavefunction is evolved through the SOFFT algorithm. 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 Figure 44. F1n at t = 0 au, for the N..II. f.i.e model. This PWTDM is evolved through the EPQCLE method. I I I 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 20 Figure 45. T/(R) at t = 14000 au, for the N..11. f.,.e model. 25 30 35 R (au) 60 65 70 75 80 P (au) 85 90 95 Figure 46. Phase space grid points at t model. 14000 au, for the N..II. f.T.e 100 0 2000 4000 6000 8000 10000 12000 14000 time (au) Figure 47. Na populations p]q and rq2 vs. time. 1 I ii SOFFT: rii / EPQCLE: 2 SOFFT: q2 0.8 "0 S0.6 .5 0.4 o 0.2 2000 nm1 .500 0 1000 0 time (au) Figure 48. Li popul:lltin r and 'l2 vs. time. 0.4 EPQCLE SOFFT 0.3  0.2 0.1 0 0.1 0.2 0.3 0.4 0 2000 4000 6000 8000 10000 12000 14000 time (au) Figure 4 9. Coherence described by Re(r/12) vs. time, for the N.. u1 f.,.:e sys tem. 0.45 uI I 0.4 0.35 0 0.3 S 0.25 S 0.2 o 0.15 0.1 0.05 0 0.05 0 500 1000 1500 2000 time (au) Figure 410. Coherence described by Re(lq12) vs. time, for the Lisurface sys tem. 0 2000 4000 6000 8000 10000 12000 14000 time (au) Figure 411. Expectation of position and dispersion for the N..i f.,.:e sys tem. EPQCLE:
