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 Wavefunction-Based Approaches ....... ..... ... 5
1.2.2 Density Operator Approaches ...... .......... 7
1.3 Quantum-Classical Liouville Equation ... . 8
1.4 Our Approach ....... ......... ... ...... 10
1.5 Simple One-Dimensional Two-State Models . ... 11
1.6 Lithium-Helium Clusters ................ .... .. 11
1.7 Outline of the Dissertation ................. .. 12
2 QUANTUM-CLASSICAL LIOUVILLE EQUATION: FORMULATION 14
2.1 Introduction .................. ........... .. 14
2.2 Quantum Liouville Equation ............ .. .. 14
2.3 W igner Representation .................. ... 15
2.4 Quantum-Classical 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 Relax-and-Drive 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 ONE-DIMENSIONAL TWO-STATE MODELS .. ...........
4.1 Introduction . . . .
4.2 Effective-Potential QCLE in the Diabatic Representation .....
4.3 Near-Resonant 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 ATOM-RARE GAS CLUSTERS: GENERAL FORMULATION 80
Introduction . . .
Physical System .. ... .. .. .. .. ..
Properties of Interest ...........
Hamiltonian for Alkali-Rare Gas Pairs ...
Hamiltonian for the Alkali-Rare 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 LITHIUM-HELIUM CLUSTERS .........
6.1 Introduction . . .
6.2 Description of the System ..........
6.3 Properties to be Investigated ........
6.4 Preparation of Lithium-Helium Clusters ..
6.4.1 Bulk Helium ... .. .. .. .. ..
6.4.2 Liquid Helium Droplets .......
6.4.3 Lithium-Helium 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 Quantum-Classical Liouville Equation 140
One-Dimensional Two-State Models ... . 141
Alkali-Rare 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 OPERATOR-FAST 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
4-1 Parameters used in the N..-- n f. .:e and Li-surface models. ...... ..42
4-2 Parameters used in the dual avoided crossing collision model. 60
4-3 Parameters used in the Nal complex model. ............ ..68
5-1 Pseudopotential rotation for d-function mixing. ..... 98
6-1 Parameters for the He-He interaction from Aziz (VA). ... 106
6-2 Parameters for the correction to the He-He interaction (V) 106
6-3 Parameters for the e--Li interaction. ................. ..113
6-4 Parameters for the e--He interaction. ................ ..113
6-5 Parameters for the Li-He core interaction. .. . ..... 114
LIST OF FIGURES
Figure page
4-1 Potential curves for Hamiltonian I: Na incident upon a metal surface. 43
4-2 Potential curves for Hamiltonian II: Li incident upon a metal surface. 44
4-3 T(R) at t = 0 au, for the N..--ii if.,.e model. This wavefunction is
evolved through the SO-FFT algorithm. ............. ..48
4-4 F11 at t = 0 au, for the N..--II. f.i.:e model. This PWTDM is evolved
through the EP-QCLE method. ................. 49
4-5 T(R) at t = 14000 au, for the N..--11. f.,.!e model. . 50
4-6 Phase space grid points at t = 14000 au, for the N..--iii f.,.e model. 51
4-7 Na populations qi and 7q2 vs. time. ................. 52
4-8 Li populations mq and 7q2 vs. time. ................. 53
4-9 Coherence described by Re(iq12) vs. time, for the ..- iii f.i:e system. 54
4-10 Coherence described by Re(M12) vs. time, for the Li-surface system. 55
4-11 Expectation of position and dispersion for the N..-- ,i f.,.:e system. .56
4-12 Expectation of momentum and dispersion for the N..--i if.i.:e system. 57
4-13 Density function p(R) for the N..-- ii f.,:e system. . ... 58
4-14 Potential curves for the dual avoided crossing collision. ... 61
4-15 Populations piq and 7q2 vs. time for the dual crossing collision model. 64
4-16 Coherence described by Re(r12) vs. time, for the dual crossing collision
m odel. ........... ...... ........ ..... 65
4-17 Grid deformation at t = 1400 au, for the dual crossing collision model. 66
4-18 Probability of transmission in the ground state, for the dual crossing
collision model. .................. ..... 67
4-19 Potential curves for the Nal complex. ................. 69
4-20 Ionic and neutral populations over time, for the Nal complex. 73
4-21 Expectation of position and its deviance, for the Nal complex. 74
4-22 Coherence as a function of time, for the Nal complex. ... 75
4-23 Phase space grid at the end of the simulation, for the Nal complex. 76
4-24 Number of steps required by the relax-and-drive algorithm, compared
to an estimated number required for a fixed timestep version. 77
6-1 Schematic of Li(2p) above a He surface. A) Li(2pw). B) Li(2p(r). 101
6-2 Radial distribution functions for bulk liquid helium .... 106
6-3 Comparison of the Aziz potential with the effective form. ...... .107
6-4 Effective He-He potential. ........ . ..... 108
6-5 Constraining potential used to keep He atoms from evaporating. 110
6-6 Temperature fluctuations of the He droplet over time .... 111
6-7 Helium density profile from the center-of-mass of the cluster. 112
6-8 Adiabatic energy for Li and He as a function of internuclear distance. 114
6-9 Adiabatic energies for Li and one or more He along the z-axis. 116
6-10 Adiabatic energy for Li and one or more He along the y-axis. 117
6-11 Adiabatic energy for Li and a surface of He atoms parallel to the x-y
plane. . . .. .. .......118
6-12 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
6-13 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
6-14 Electronic population of Li as it emerges from the He cluster. 123
6-15 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
6-16 Mixing of the Li(2pa) and Li(2pr) states at distances where Li(2p) is
triply degenerate. ................ ........ 126
6-17 Electronic population of Li with and without a perturbing electro-
magnetic field, resonant to the D line. . 128
6-18 Dipole emission spectra of Li(2pa) as it recedes from the He cluster
surface. .... .. .. .... ..... ...... 129
6-19 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
6-20 Electronic population of Li(2p7) as it interacts with the He cluster
surface. . .. .. .. .. 132
6-21 Dipole emission spectrum of Li(2pr) during the first 3000 au. .... 134
6-22 Dipole emission spectrum of Li(2pr) during the final 3000 au. .... 135
6-23 Adiabatic curves of Li surrounded by a cubic lattice of He atoms. The
parameter R refers to the half-length of the lattice edge. 136
6-24 Decay of Li(2pr) surrounded by surface He atoms, induced by an EM
field with frequency resonant to the Li(2pr <-- 2sra) transition. .138
A-1 Flowchart describing the cauldron program. . 148
C-1 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 quantum-classical 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 quantum-classical method is applied to the study of three simple
one-dimensional two-state models. The first model represents the photoinduced des-
orption of an alkali atom from a metal surface, where near-resonant 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 long-range 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 atom-rare gas interactions with i-dependent
semi-local 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 many-atom 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 many-body molecular systems. In order to study these systems,
we are particularly interested in mixed quantum-classical 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 quantum-classical 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 few-atom 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 alkali-rare 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 quantum-classical 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 time-invariant
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.in-i. .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 zero-point 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 Wavefunction-Based Approaches
Wavefunction-based 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 self-consistent 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 self-consistent 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 quantum-classical
schemes retain this mean-field 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 quantum-classical 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 Maxwell-Schrodinger 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 spin-boson 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, -1i--i 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 quantum-classical 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 self-consistent 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 Hellmann-Feynman 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
operator-based 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 quantum-classical schemes. In Section 1.3, we outline some of the major
approaches found in the literature.
1.3 Quantum-Classical 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. So-called 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 quantum-classical 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 quantum-classical 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 off-diagonal 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 off-diagonal 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
quantum-classical Liouville equation (EP-QCLE) 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 higher-order approximations and quantum-classical
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 trajectory-based numerical solutions that reduce to the classi-
cal Liouville equation in the classical limit, providing an intuitive computational
framework.
Our EP-QCLE 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 Hellmann-Feynman 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 EP-QCLE 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 one-dimensional models expanded in
a two-state diabetic basis, to fully three-dimensional 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 One-Dimensional Two-State Models
We test our methods on a set of one-dimensional two-state 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 near-resonant 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
long-range attraction of the excited state results in oscillatory nuclear motion.86
1.6 Lithium-Helium 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. Semi-local i-dependent
pseudopotentials are used to describe the lithium-helium interactions, because they
are found to accurately reproduce adiabatic energy curves for the lithium-helium
pair for s, p and d symmetries.90'91'92'93'94 By introducing the pseudopotential for-
malism into the EP-QCLE, 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 quantum-classical systems. The ideas are presented for the general case of
quantal and quasiclassical variables, and it is shown how the EP-QCLE naturally
emerges from the partial Wigner transform over the quasiclassical variables.
Chapter 3 explores the computational aspects of the EP-QCLE. 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 one-dimensional two-state
-,-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 alkali-rare gas clusters. We describe
the Hamiltonian for general alkali-rare 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 lithium-helium 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 EP-QCLE for the test model systems and the full
lithium-helium 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 quantum-classical dynamics with the
full quantal treatment.
CHAPTER 2
QUANTUM-CLASSICAL 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 quantum-classical meth-
ods, because the classical limit of the DOp is the classical density function. One
particular route to a mixed quantum-classical 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 classical-like (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 {l|i)} 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 quantum-classical 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/2|r + 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 quantum-classical 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 N-body
-v-1, 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' = r-z/2, we retrieve the quantum expectation
r, p)Aw(r,p) = dqdq'F(q, q')A(q', q)
= dqdq'(q|Pfq')(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 Quantum-Classical 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 electronic-nuclear
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 quantum-classical 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 (aRVa-Rv)
to zero,
Trqu 1W ( RV a VI
OR OR
(2.26)
In this way, we retrieve the Hellmann-Feynman 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 effective-potential QCLE (EP-QCLE),
0f1w 1 P 80w 0 aw
Ot i[H F M OR F P
The approximations used to derive the EP-QCLE substantially reduce the
quantal character of the quasiclassical solution space. In contrast, the best adia-
batic methods, based on the Born-Oppenheimer separation of nuclear and electronic
motion, retain full quantum nuclear dynamics along adiabatic curves.102 A princi-
pal advantage to the EP-QCLE is its nonadiabatic character, which is capable of
describing nonadiabatic events when the Born-Oppenheimer 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 EP-QCLE are discussed in
the next chapter.
CHAPTER 3
QUANTUM CLASSICAL LIOUVILLE EQUATION: COMPUTATIONAL
ASPECTS
3.1 Introduction
The EP-QCLE 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
EP-QCLE 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 EP-QCLE 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) S-l (FIIl4)S-1 (4,, (3.5)
F
A 4) S-1 (4 AI4) S- I(4 (3.6)
A
where S is the overlap (4|4). Projecting Eq. 3.3 on this basis set, and setting
h 1, we obtain
dT
r(iH q- _t)S-1 S-l(iH' + )F (3.7)
dt
where we have used the notation,
dR dP
S (fl (| //&R|4) + (4d|/&P|4). (3.8)
dt dt
3.4 Nuclear Phase Space Grid
Although we've transformed the EP-QCLE 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 Relax-and-Drive 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-
and-drive 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 relax-and-drive procedure
has been shown to give excellent results for a wide variety of two-body 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 EP-QCLE for
a single trajectory,
S =. W(t)F(t)- F(t)Wt(t), (3.15)
dt
where we have defined
W S-1(-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(t-to)]. (3.18)
If we diagonalize Wo,
Wo = TAT-1, (3.19)
we can rewrite Eq. 3.17 as,
F(t) T[T-Uo(t, to)T]T- o(Tt)-1[TtUt(t, to)(Tt)-1]Tt. (3.20)
The exponential matrices can be formed analytically, since
T-1Uot, to)T exp[-iA(t-to)], (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 relax-and-drive 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 relax-and-drive algorithm,
although to keep accuracy to O(At2) it is necessary to integrate using an algorithm
like velocity Verlet or Runge-Kutta. We choose to use the velocity Verlet method,
which is accurate to O(At2) and is self-starting. 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 relax-and-drive 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)]T-lro(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 = T-1D(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)T-lro(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)S-1/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=)S-1/2r'S-1/2(4
) | F)r 1|. (3.53)
Thus we equate
S S-1/2F'S-1/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)S-1/2A'S-1/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[S-1/2AwS-1/2S1/2FrS1/2]
STrqu[AwFw]. (3.60)
3.6.4 Hamiltonian Eigenstates and Eigenvalues
Using the orthonormal representation of the Hamiltonian,
H' = S-1/2HS-1/2, (3.61)
we can compute its eigenvalues by diagonalizing the matrix,
L-1H'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 object-oriented 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 object-oriented 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 relax-and-drive, 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
,-i-1-. 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
ONE-DIMENSIONAL TWO-STATE MODELS
4.1 Introduction
While our ultimate goal is to study realistic three-dimensional 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 near-resonant 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 EP-QCLE approach to the dynamics of the full
quantal system. For all three models, we will see that the mixed quantum-classical
methods provide qualitatively, and often quantitatively similar results to the full
quantal evolution.
4.2 Effective-Potential QCLE in the Diabatic Representation
In Chapter 3, we derived the EP-QCLE 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
EP-QCLE 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 (SO-FFT) 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 SO-FFT procedure provides an excellent test of the accu-
racy of the EP-QCLE. A complete description of the SO-FFT is given in Appendix
B, and the code implementing the SO-FFT (qualdron) is described in Appendix C.
4.3 Near-Resonant 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 near-resonant
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 4-1. Parameters used in the N..--U1 f...:e and Li-surface 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 Ro-R 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 4-1, where we con-
sider two model Hamiltonians: (I) N..--II[ f.,:e and (II) Li-surface. The diabetic
potentials for Hamiltonians I and II are shown in Figures 4-1 and 4-2.
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 4-1. 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 4-2. 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 off-diagonal 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 quantum-classical 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
SO-FFT 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 EP-QCLE
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 4-3 and 4-4 show the initial wavefunction and its PWT, respectively.
Note that the PWTDM formed from a Gaussian wavepacket is a Gaussian function
itself, albeit in two-dimensional phase space. Figure 4-5 shows the wavefunction at
t = 14000 au, having been propagated through the SO-FFT method. The PWTDM
evolves in phase space through the EP-QCLE, and in Figure 4-6 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 .i-i. ,,,.1, -. i' state, where each point is subject to a
vanishing Hellmann-Feynman force, and thus propagates at constant velocity.
The observables are presented in Figures 4-7 to 4-13. There are a number of
interesting things we can glean from these plots. From Figures 4-7 and 4-8 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 4-9 and 4-10 describe the coherence between the states. For the
N.- -ii f.i.e system, the coherence remains large for long times, while for the Li-surface
3 3.5 4 4.5 5 5.5 6 6.5 7
R (au)
Figure 4-3. T(R) at t = 0 au, for the N..--i i f.-.e model. This wavefunction
is evolved through the SO-FFT algorithm.
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
Figure 4-4. F1n at t = 0 au, for the N..--II. f.i.e model. This PWTDM is
evolved through the EP-QCLE method.
I I I
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
20
Figure 4-5. 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 4-6.
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 4-7. Na populations p]q and rq2 vs. time.
1 I ii SO-FFT: rii
/ EP-QCLE: 2
SO-FFT: q2
0.8
"0
S0.6
.5 0.4
o
0.2
2000
n-m1 .500
0 1000
0 time (au)
Figure 4-8. Li popul:lltin r and 'l2 vs. time.
0.4
EP-QCLE
SO-FFT
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 u-I- 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 4-10. Coherence described by Re(lq12) vs. time, for the Li-surface sys-
tem.
0 2000
4000
6000 8000 10000 12000 14000
time (au)
Figure 4-11. Expectation of position and dispersion for the N..--i f.,.:e sys-
tem.
EP-QCLE:
SO-FFT:
EP-QCLE: Gp
SO-FFT: op
0 2000 4000
6000
8000 10000 12000 14000
time (au)
Figure 4-12. Expectation of momentum and dispersion for the N..--i fi .e
system.
58
0.16
EP-QCLE
SO-FFT
0.14
0.12
0.1
- 0.08
0.06
0.04 -
0.02 -
0
0 5 10 15 20 25 30 35 40
R (au)
Figure 4-13. Density function p(R) for the N..--U1 f[i.:e system.
-i,-I-1 iii. the coherence rapidly diminishes once the interaction potential is crossed.
This can be explained by the size of the energy gap at ..-1i.. .1l i' distances, as
larger ionization energies lead to more rapidly vanishing coherences.
Figure 4-11 shows the expectation of the position of Na steadily increases
from its initial average, while the dispersion in position initially decreases and then
increases again. On the other hand, Figure 4-12 shows a marked difference in the
behavior of the expectation of the momentum, where it begins at (P) = 0 au, rapidly
increases and then becomes stationary around (P) = 78 au. We also find that the
momentum dispersion first increases, and then decreases.
Finally, Figure 4-13 presents the density function at the end of the simulation.
Because of the distortion of the phase space grid, the density function obtained
from the EP-QCLE had to be calculated using bins and through approximations
within these bins, and thus is subject to some noise. Nevertheless, we find excellent
agreement between the EP-QCLE and SO-FFT results.
For all observables, the EP-QCLE results are quantitatively similar to the SO-
FFT results to visual resolution. Having studied both Na and Li approaching a
metal surface, we see that the EP-QCLE can be expected to yield very accurate
results for these kinds of systems, even when coherence is maintained over long
periods.
4.4 Binary Collision Involving Two Avoided Crossings
4.4.1 Model Details
For the second system, we look at a two-state collision model where the diabetic
surfaces intersect twice. Because of the dual crossing and the coupling in this region,
quantum interference and effects such as tunnelling play a substantial role in the
dynamics of the quasiclassical variable. As such, this model is quite demanding
for mixed quantum-classical methods, where one can expect deviations at lower
energies.
Table 4-2. Parameters used in the dual avoided crossing collision model.
Parameter Value (au)
Ro -8
Po [10, 30]
a 2.5176
M 2000
The Hamiltonian elements are79
H1 (R) = 0, (4.22)
H22(R) = -0.1exp(-0.28R2) + 0.05, (4.23)
H12(R) = 0.015 exp(-0.06R2). (4.24)
In the above, R is the nuclear-nuclear separation. The initial state is a ground state
Gaussian wavepacket which begins in the .i-.,mptotic region R = 8 au. Its PWT
over R gives the Gaussian density,
ll(P, R) = 1 [ R) 2p 2 (4.25)
with F12 F21 F22 = 0. At t 0, the wavepacket propagates toward the
region of coupling, where it is partially transmitted and partially reflected, now
with populations in both the ground and excited state. The parameters used in the
calculation are shown in Table 4-2. The diabetic potentials used in the model are
shown in Figure 4-14.
4.4.2 Properties of Interest
As in the alkali atom-surface model, we can follow the populations over time.
The system begins in the ground state, and as it passes through the collision region,
we expect some of the population to transfer to the excited state. The amount of
transfer depends on the collision energy. We can also follow the coherence once the
H -
H22
H12
-
-10 -5 0 5 10
R (au)
Figure 4-14. Potential curves for the dual avoided crossing collision.
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
collision has ended and the system reaches ..i mptotic values. We also consider the
transmission of the ground state as the wavepacket passes through the interaction.
Transmission. As the wavepacket passes through the collision, part of it
continues forward while the remainder is reflected back. We can study this by
computing the probability of transmission of the ground state,
0
^ = dRdPTn{RP). (4.26)
Although not included, we could also compute the probability of reflection of the
ground state, as well as the .1- i.l.1 .tic populations of the excited state.
4.4.3 Results
In Figure 4-15 we see the total population transfer at energy Po = 30 au. In
this energy region, the EP-QCLE matches the quantum results to within a cou-
ple percent. For both the EP-QCLE and SO-FFT simulations, the reflection of
the wavefunction (or PWTDM) was negligible, so that we need only compare the
transmission. This negligible reflection was found throughout the energies studied.
In Figure 4-16 we see the coherence as a function of time. The variations are
large initially and diminish .i ,'iiil' .tically to regular oscillations. The EP-QCLE
captures this behavior qualitatively, and to within 5% quantitatively.
We show the deformation of the phase space grid in Figure 4-17. We see
that it is primarily the inner part of the grid that undergoes deformation, while
the enclosing points maintain a structured order. This is likely due to only a small
range of values in phase space which are seriously affected by the potentials. Outside
the interaction region, all Hellmann-Feynman forces are zero, so one would expect
deformations only for points whose time in the interaction region was significant.
In Figure 4-18 we display the transmission probability for a wide range of ener-
gies. As expected, the EP-QCLE performs better for higher energies, and while the
double-well is qualitatively reproduced, the EP-QCLE fails to give quantitatively
accurate results for energies lower than (log E = -2). This is likely due to the quan-
tum tunnelling effects described earlier, which are not expected to be reproduced
well by the EP-QCLE. We also compare transmission probabilities obtained through
surface hopping methods by Tully and coworkers.51 These surface hopping calcu-
lations show deviations from the quantal results that are similar to the EP-QCLE
probabilities at lower energies; at higher energies, the EP-QCLE is slightly superior
to the surface hopping scheme.
4.5 Photoinduced Dissociation of a Diatomic System
4.5.1 Model Details
For the third test system, we explore a model of the Nal complex. As in
the previous models, it involves two diabetic surfaces and an interaction around
the avoided crossing. A substantially different feature, however, is a long-range
Coulombic attraction in an ionic state. As we shall see, this attractive potential
results in the the complex oscillating between a neutral and ionic state as the sodium
and iodine separate and come back together, partially dissociating at each crossing
into an .i-, ,iii'l.itically neutral state. This oscillatory motion is a good test for the
EP-QCLE at long times in cases where ..i- mptotic states are not reached quickly.
The Hamiltonian elements .i "''
H11(R) = Aiexp[-/3(R- Ro)], (4.27)
H22R) [A2 (B2/R)] exp(-R/p) /R- (+ -)/2R4
-C2/R6 2+-/R7 + AEo, (4.28)
H12(R) A12 exp[- 12(R R)2]. (4.29)
To form the initial state, we Taylor expand (to first order) the ionic potential H22
about its minimum R = Ro, and find the lowest energy state of this harmonic well.
At t = 0, the wavepacket undergoes a sudden optical promotion to the neutral curve,
EP-QCLE: rl
SO-FFT: l
EP-QCLE: 12
SO-FFT: 12
0 200 400 600 800 1000 1200 1400
time (au)
Figure 4-15. Populations qrl and rl2 vs. time for the dual crossing collision
model.
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2
0 200 400
Figure 4-16.
EP-QCLE
SO-FFT
600 800 1000 1200 1400
time (au)
Coherence described by Re(riq2) vs. time, for the dual crossing
collision model.
ue 4-1
(rid ,, rMu/ 3) 34
S400 a, for le dualoss
0coli,..i,
67
S" / \ EP-QCLE
S------- SO-FFT -
S0.9- Tully et al. -
0.8
0.7
I- \\
F-
0.6 -
0.5
0.4
0.3
-4 -3.5 -3 -2.5 -2 -1.5 -1
loge(E) (au)
Figure 4-18. Probability of transmission in the ground state, for the dual
crossing collision model.
Table 4-3. Parameters used in the Nal complex model.
Parameter Value (au)
Ro 5.047
Po 0.0
a 0.12462
A1 0 i111 i ',
A2 101.43
A12 0.00202
B2 3.000
C2 18.950
A+ 2.756
A- 12.179
p 0.660
AEo 0.07626
/1 2.158
/12 0.194
R, 13.24
M 35,482
so that the PWTDM becomes,
rl(P, R) 1 R RO 2 2(p o)2 (4.30)
with F12 F21 F22 = 0. The simulation follows the resulting motion of this
state. The model parameters are shown in Table 4-3. The diabetic potentials are
displayed in Figure 4-19.
4.5.2 Properties of Interest
In addition to observing the coherence, expectation value of the position and its
dispersion, and the phase space grid, we also consider bound and free neutral and
ionic populations as the Nal complex oscillates from its primarily ionic to primarily
covalent state.
Bound and free populations. It is interesting to follow the populations of
the ionic and neutral states as the system evolves, giving insight into the nature of
the dissociation into an ..I-.i.l.1 i ..lly free neutral system. For this, we define three
5 10 15 20 25
R (au)
Figure 4-19. Potential curves for the Nal complex.
0.05
0
-0.05
populations: the ionic, the bound neutral and the free neutral. The ionic popula-
tion is the probability of finding the system in the ionic state at any internuclear
separation,
00
/2 = dRdP22(R, P). (4.31)
-00
We define the bound neutral population is the probability of finding the Nal complex
in the neutral state at nuclear separation up to the crossing point R = R,,
R
f = fdRdPI (R, P), (4.32)
0
while the free ionic population is defined as the probability of the neutral state from
the crossing point beyond,
J{ = dRdPFi(R, P). (4.33)
R
This division between bound and free is motivated by the observation that the
majority of the transfer between ionic and neutral potentials occurs at the avoided
crossing, and that any part of the wavepacket which ends up in the neutral state
but propagating toward infinite separation has a negligible probability of shifting to
the ionic state much beyond the avoided crossing.
4.5.3 Results
The ionic and covalent populations are displayed in Figure 4-20. We see the
oscillations in the populations between ionic and covalent, repeating approximately
every 40 000 au. This pattern can be compared to the expectation values of position
see in Figure 4-21. The position oscillates with the same frequency as the change
in population, showing that each time the wavepacket heads across the avoided
crossing from its bound covalent state, it converts almost completely into the ionic
state, with a small amount escaping into the free neutral state. Over time, the
free neutral population gradually increases at these crossings, and the Nal slowly
dissociates. The EP-QCLE is quantitatively similar to the exact results for the first
half of the simulation, and maintains qualitative accuracy for the remainder.
The dispersion in Figure 4-21 shows an interesting difference between the exact
and the quantum-classical algorithm. While the dispersion continues to rise in the
SO-FFT simulation, it reaches its first peak and then begins to decline somewhat in
the EP-QCLE simulation. This reflects the nature of the effective potential, where
each point is guided by a combination of excited and ground state forces. Because
the ionic curve does not permit escape, what would normally be ..i- .mptotically free
wavepackets tend to be pulled back toward the crossing because of the attractive
ionic potential. Consequently, the free neutral population is always lower using the
EP-QCLE equation than the SO-FFT, an observation supported by Figure 4-20.
Since the majority of the population remains in the ionic or bound neutral state,
and these populations are well matched between the exact and quantum-classical
simulations, it is not surprising that the expectation value of the position is quan-
titatively in the beginning, and qualitatively for the i. is,.iinii.-. similar for both
methods. However, the dispersion is much more sensitive to the increased ..I-i.mp-
totically free wavepackets in the SO-FFT simulation, and for the reasons discussed,
we find significant divergence between the EP-QCLE and SO-FFT results.
The coherence, shown in Figure 4-22 initially peaks through the first crossing,
but through subsequent crossings it is substantially diminished. However, the EP-
QCLE shows quantitatively similar results to the SO-FFT calculations.
The deformation of the phase space grid, plotted in Figure 4-23, has char-
acteristics not seen in the other two models. One line of points emits from the
center of the cluster, quickly straightening and reflecting the negligible force on the
points. These points correspond to the .,i- ,,iiilitically free neutral components of
the PWTDM. The second group circles around, gaining velocity and position, then
turning. These ellipses are characteristic of the phase space of classical particles in
a well, and indeed reflect the quasiclassical motion under the Hellmann-Feynman
force of the PWTDM points as they follow the ionic and bound covalent curves.
4.6 Comparison Using Variable and Constant Timesteps
In this Section, we evaluate the usefulness of the variable timestep aspect of
the relax-and-drive algorithm. To do this, we consider the N..--II f.,:e algorithm of
Section 4.3, and simulate using varying upper and lower tolerances. The number
of steps taken in each case is compared with the number that would be required of
the same algorithm, but keeping the timestep fixed. The fixed timestep would nec-
essarily advance by steps no greater than the smallest timestep used in the variable
timestep algorithm, and it is based on this timestep that we estimate the corre-
sponding steps required for the fixed timestep approach.
The results are shown in Figure 4-24. We see that as the tolerance decreases
(and thus the accuracy increases), the fraction of steps saved by the introduction of
the variable timestep increases superlinearly. One concludes that while the variable
timestep may not be important for low accuracy simulations, considerable compu-
tational savings can be had for high accuracy propagation.
4.7 Conclusion
By examining three simple two-state models, we were able to compare the
EP-QCLE method with the exact SO-FFT quantum mechanical solution. For all
models, we found very good agreement between the EP-QCLE and the SO-FFT
results. The agreement was at least qualitative, and in many cases quantitative to
visual precision. We also saw conditions under which the quantum-classical model
deviated from exact quantal results. Finally, we compared a fixed timestep variant
of the relax-and-drive algorithm with the variable timestep version.
For the model of an alkali atom approaching a metal surface, we examined prob-
ability transitions, position expectation and its deviation, momentum expectation
73
2
EP-QCLE: Ionic
SO-FFT: Ionic
EP-QCLE: Neutral Bound
SO-FFT: Neutral Bound
EP-QCLE: Neutral Free
C 1.5 SO-FFT: Neutral Free
0
0
0
0z
.--" ------
"oz- -- -- -- --
o 0.5 .
-- --- --- -
0 T-- --- T- ---------- I-- I'^ 'T '''
0 20000 40000 60000 80000 100000 120000
time (au)
Figure 4-20. Ionic and neutral populations over time, for the Nal complex.
74
25
EP-QCLE:
SO-FFT: -:i-
EP-QCLE: GR
SO-EFT: o -:
20
15
A
v
10 / -.---"
5
0 -I I I I I I I I
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 100000
time (au)
Figure 4-21. Expectation of position and its deviance, for the Nal complex.
0.35
EP-QCLE
SO-FFT
0.3
0.25
(1)
0 0.2
rr
o
| 0.15
0
t-
0.1
0.05 -
0
0 20000 40000 60000 80000 100000 120000
time (au)
Figure 4-22. Coherence as a function of time, for the Nal complex.
140
120
100
80
rr
60
40
20 ..
0
-80 -60 -40 -20 0 20 40 60 80
P (au)
Figure 4-23. Phase space grid at the end of the simulation, for the Nal
complex.
3500
Variable
Fixed --x
3000
2500
0-
0)
a 2000
E
0
5 1500 -
E
1000
500 -
103-10-1 10-5-10-3 10 -5 10-9- 10-7
tolerance
Figure 4-24. Number of steps required by the relax-and-drive algorithm,
compared to an estimated number required for a fixed timestep
version.
and its deviation, density function and phase space grid evolution. We found that
the EP-QCLE method reproduced quantitatively the exact values found through the
quantum model for all observables. We also found the phase space grid to deform
substantially as the system evolved, initially distorting but eventually lining up in
a straight line as the effective force vanished.
For the dual crossing diabetic surface collision, in addition to probability tran-
sitions for a given energy, we calculated the transmission probability of the ground
state for a wide range of energies. This showed the EP-QCLE to deviate at the lower
energies, where nuclear interference effects were important. At higher energies, the
correspondence between the mixed quantum-classical results and the exact quantum
results was good. We also noted that the phase space grid deformed only for points
which spent a significant amount of time in the region of strong interaction. This is
reasonable, as the deformation comes from the effective force, which is non-varying
(and in fact zero) outside the interaction.
The Nal model showed interesting back-and-forth transitions between the ionic
and covalent state, as the internuclear distance oscillated due to the long-range
attractive ionic potential. We also saw that each passing through the avoided cross-
ing led to a small amount of dissociation into the ..i-1.mptotically free neutral state.
Ionic, bound neutral and free neutral populations were quantitatively similar for the
EP-QCLE and SO-FFT algorithms for the beginning of the simulations, and quali-
tatively similar for the remainder. Coherences were quantitatively similar, showing
an initial peak after the first pass through the avoided crossing, and much smaller
magnitudes thereafter. Finally, the grid deformation showed the grid to split into
two groups. The first represented the .i -1'mptotically free neutral state, where the
points form a straight line as the Hellmann-Feynman force vanishes and the parti-
cles propagate to infinite distances. The second group circled in phase space, as the
79
points oscillated back and forth in the potential well formed by the ionic attraction
at large distances and the neutral repulsion at small distances.
The EP-QCLE method was shown to be very robust under a wide variety of
conditions. Comparing its variable timestep version (which was used in all the simu-
lations) to the fixed timestep alternative, we found that not only was the EP-QCLE
method accurate, but use of the relax-and-drive propagation generated much more
efficient simulations than could have been obtained without dynamically varying the
timestep.
CHAPTER 5
ALKALI ATOM-RARE GAS CLUSTERS: GENERAL FORMULATION
5.1 Introduction
In this chapter, we describe the interactions of an alkali atom (Ak) embedded in
a cluster of rare gas atoms (RgN). Cluster dynamics are interesting insofar as they
present a bridge between isolated atoms and the (effectively) infinite bulk liquid, and
can provide insight into the dynamics of molecules on extended surfaces. The Ak-
RgN clusters are among the simplest cluster systems to study because of the presence
of a single valence electron in the alkali atom, and the closed shell structure of the
rare gas atoms. Pseudopotential descriptions of the interactions between the valence
electron (e-), the Ak core and the Rg atoms can greatly simplify the description of
the cluster with very modest penalties in accuracy. In the context of mixed quantum-
classical models, these clusters are well suited to simulation through the EP-QCLE
by treating the valence electron quantally and the nuclear cores quasiclassically. For
the majority of the chapter, we consider the general case of an alkali atom embedded
in a rare gas cluster, although our simulations in Chapter 6 focus on the specific
case of the lithium atom (Li) embedded in a helium cluster (HeN).
5.2 Physical System
We consider a cluster initially at thermal equilibrium, and introduce a ground
or excited alkali atom to the center of the cluster. We do not concern ourselves
with the means by which the alkali atom is excited or embedded within the cluster,
although typically the excitation is due to a laser pulse.
From this initial setup, we follow the dynamics of the alkali atom as it deexcites
and moves within the cluster. We further assume the cluster is isolated from any
environmental effects, so that spectral and configuration measurements represent
those of an isolated Ak-RgN cluster. As we see later, we introduce a containing
potential to keep the Rg atoms from dispersing. This potential is used strictly to
maintain a cluster formation, and does not represent interaction with an environ-
ment.
5.3 Properties of Interest
First of all, we are interested in the structure of the Rg cluster in equilibrium,
before the introduction of the Ak atom. We need to reproduce known density profiles
and pair-pair correlation functions to ensure the cluster is representative of a real
ll,-i, .1i system. Secondly, we are interested in following the dynamics of the Ak
atom as it moves within the cluster. In particular, we are looking for migration from
the center of the cluster to its surface, as other experimental and theoretical studies
indicate that the Ak atom tends to reside on the cluster surface. Finally, we wish
to compute the time dependent emission spectra resulting from the deexcitation of
the Ak atom,
Ak(n'l') +RgN Ak(nl) + RgN + (5.1)
where 0 is a photon with energy corresponding to the electronic decay.
5.4 Hamiltonian for Alkali-Rare Gas Pairs
Following the dynamics of an Ak-RgN cluster through a full quantal treatment
would be an extremely demanding problem computationally for more than a few Rg
atoms. Instead, we reduce the Ak-Rg interaction to a three-body problem by using
pseudopotential interactions between the nuclear cores and the electron. The pseu-
dopotential treatment is explored extensively for the single Ak-Rg pair in Reyes.109
In this section, we summarize this approach.
In order to describe the three-body interaction, we consider a fixed nuclear
configuration with RAB the position vector from the alkali core (A) to the rare gas
atom (B), and rA (rB) the position vector from the alkali core (rare gas atom) to
the electron e-. With this notation, we can write the Hamiltonian in five distinct
components,
pair 2 r Ak(r)
f- = -VA A A)+VBgrB)+VA cro A, AB)
+V1ore(RAB). (5.2)
We have implicitly taken the PWT over nuclear variables but not electronic variables.
Because the potentials can all be expressed as polynomial functions of Q, the PWT
amounts to replacing Q by R throughout the Hamiltonian. In what follows, we
will assume we are working with partially Wigner transformed operators, but drop
the subscript 'W' for notational simplicity. Also note that the electronic variable
remains quantal, although we will use the notation r rather than q for consistency
with common usage in the literature.
The first term on the RHS of Eq. 5.2 is, of course, the kinetic energy operator
of the valence electron. The second and third terms are the potentials arising from
the interaction of the electron with the Ak core and the Rg atom, respectively. The
VA'O term is a crossterm stemming from the polarization of the Rg atom by both
the valence electron and Ak core. Finally, the last term Vo" is the interaction
between the Ak core and the Rg atom. We examine each of these potentials in
detail for the general Ak-Rg pair, and provide specific parameters for the Li-He
interaction in Chapter 6.
The e--Ak core potential can be divided into three components,
V k (AZA ^tAk
k(rAk ) -Z VA (-A) V+T(rA), (5.3)
TA
where the first term is the Coulomb interaction between the Ak core of charge ZA
and the electron. The second contribution arises from the dipole polarization of the
Ak core by the valence electron,
d
VT^'j -A A) A- w(A, 6A)2, (5.4)
2rA
with a$ the dipole polarizability of the core and w a cutoff function of the distance
6,
w(r, 6) = exp(-6r2). (5.5)
The final contribution to VAk is an i-dependent short-range pseudopotential,
r (rx) B1, exp(- / )P1x, (5.6)
l,i
where B,,i and /3,i are pseudopotential parameters adjusted to fit experimental data,
Pix is the projection operator on angular symmetry 1,
Pix = Y Yim(rx))(Yim(x)I, (5.7)
and the states { Yim(x))} are the spherical harmonic functions centered on the core
X. This pseudopotential simulates repulsion due to the effects of the Pauli exclusion
principle when the valence electron approaches the core electrons, as well as the
attraction due to the incomplete screening of the nuclear charge.
The third term in Eq. 5.2 is the e--Rg potential, and can be divided into two
terms,
V""(rB) pol-(r) + *(rB). (5.8)
v g^rB VToTM (5.8)
The first component stems from the polarization of the Rg atom by the electron,
d q4
tR(r) (rB, 6)4 (r, 6) (5.9)
VB ()r
where aC a/ 6/1, a' being the quadrupole polarizability of the Rg atom
and -6/1 the dynamical correction to the static polarizability. The short-range
pseudopotential VJ' is defined in Eq. 5.6.
The fourth term in Eq. 5.2 is a crossterm arising from the polarization of the
Rg atom by both the Ak core and the valence electron,
Ir Pi (cos) q P2 PCos0)
VAo(rA, RAB) ZA 2 2 WB, 6)2 ZAa 2 2 -B 6)3. (5.10)
IAB B 1)ABrB
Here, P1 and P2 are the Legendre polynomials and 0 is the angle between the vectors
RAB and rB.
Finally, the last term in Eq. 5.2 is the interaction between the alkali core and
the rare gas atom. It is assumed to involve short-range, dipole and quadrupole
contributions,
I1 1 '1
VA re(RAB) r V(RAB) 2 B 2 2 BB 3 (5.11)
where ac' = a + 2a' d2. The short-range term VjI' has been determined by
assuming it to have the form,
VI (RAB) = Aexp(-bRAB), (5.12)
and adjusting the parameters A and b so that VAI'e(RAB) in Eq. 5.11 fits the most
accurate curves in the literature. This procedure has been shown to accurately
reproduce Ak-Rg energy curves for a variety of alkali atom and rare gas combina-
tions.109
5.5 Hamiltonian for the Alkali-Rare Gas Cluster
The cluster is similar to the pair, in that each Ak-Rg pair in the cluster is
treated as a three-body problem involving alkali core, valence electron and rare
gas atom, and described using the same Hamiltonian elements. In the case of the
cluster, however, we also require that the Rg atoms maintain a cohesive structure
rather than drift apart. To ensure this, we impose a constraining potential on the
;- 1 I III
Defining RABi as the position vector from the alkali core to the Rg atom i, and
RB, the position vector of the Rg atom from the origin, we can write the cluster
Hamiltonian as,
icluster 1 2 Ak NR B
i=1
2 A (rA) + vY, rB)
N
i=1
N
+ VCore(R RB, ). (5.13)
j>i=1
The final term in Eq. 5.13 is the interaction between the Rg atoms for the
quasiclassical motion. Following Aziz and coworkers,110 we take the Rg-Rg potential
to have the form, for internuclear distance R,
Veore(R) = V*(x), (5.14)
where
V*(x) A* exp(-a*x + *2) F(x) 2 2 (5.15)
2 -;+ '6
j=0
F (x) = exp[-(D/x 1)2], x < D (5.16)
F(x) (5.16)
1, x >D
This potential has been written in terms of the dimensionless distance x = R/RM.
The constraining potential, Vhold is taken to be a sigmoidal function centered
along the boundary of a sphere of radius Rhold,
hold (RB ) a ex(5.17)
1+ exp [-b(RB, Rhold)]'
This potential shows .'.-mptotes Vhold(-oo) = 0 and Vzold(+oo) = a, with midpoint
at the holding radius Vhod(Rhold) = a/2. The steepness of the function is controlled
by the parameter b, and determines the strength and range of the holding force. It
acts on the Rg atoms only, keeping them bound roughly within a sphere of radius
Rhold, while permitting the alkali atom free motion within the cluster.
5.6 Electronic Spectral Calculations
When a molecular system undergoes electronic motion, the accelerating charges
emit electromagnetic radiation. At distances large compared to the electronic motion,
the flux of this radiation at point r from the center of the system is given by the
Poynting vector,111
1
S(t) -(Ex B)
Pto
to 2sin2 2
I-i[t)2 r (5.18)
where D(t) is the dipole moment of the system, [to is the permittivity constant and
c is the speed of light. By integrating over all angles, we obtain the power emited
by the source,
P(t) S(t) da
I |A (t)|2 (5.19)
where
A(t) (t). (5.20)
From Eq. 5.20, we can see that the computation of the dipole moment is essential
to the spectral calculation. In a mixed quantum-classical system, the dipole moment
is obtained by calculating the expectation value of the dipole operator D,
D(t) Tr[F(t)D]/Tr[P(t)].
(5.21)
