Information • Textbooks • Media • Resources edited by
teaching with technology
James P. Birk Arizona State University Tempe, AZ 85287
A Spreadsheet Template for Quantum Mechanical Wavepacket Propagation John C. Hansen Chemistry Department, Southwest State University, Marshall, MN 56258 Donald J. Kouri Departments of Chemistry and Physics, University of Houston, Houston, TX 77204-5641 David K. Hoffman Department of Chemistry and Ames Laboratory, Iowa State University, Ames, IA 50011 In the typical physical chemistry course, one semester is devoted to quantum mechanics and its chemical applications. Usually these courses deal exclusively with solutions to the time-independent Schrödinger equation for systems of chemical interest. If the time-dependent Schrödinger equation is mentioned at all, it is usually for the purpose of carrying out a formal separation of variables to derive the eigenvalue equation for the energy levels. Although the determination of bound-state energy levels certainly covers the most common chemical applications, a number of problems of chemical interest must be omitted if time-dependent phenomena are completely ignored. For example, chemical reaction dynamics are most naturally described using the time evolution of the system, and they typically require significantly more sophisticated mathematical techniques when discussed from a time-independent quantum mechanical point of view. Furthermore, the exclusive concentration on stationary states in physical chemistry may leave students with serious misconceptions about quantum mechanics. Although chemistry students are familiar and comfortable with the connection between wave functions and probability distributions, they are not familiar with wave functions as dynamically evolving entities describing a particle’s motion. Most are unaware of the intimate connection between quantum and classical mechanics, which is most apparent in the analysis of quantum mechanical wavepacket motion. Some may not even be aware that time-dependent wave functions exist. The usual emphases in a physical chemistry course can easily leave students with the impression that stationary states are the only possible states. Hence they would have no sense of the dynamical implications of quantum mechanics. Although it may be desirable to give students an exposure to time-dependent phenomena in physical chemistry, one semester is simply not enough time for any kind of indepth theoretical discussion. A practical alternative may be to present several numerical solutions in order to show the qualitative behavior of time-dependent wave functions. A number of topics in an introductory physical chemistry course can be understood or illustrated using simple, timedependent quantum scattering. Examples include reflection/transmission in scattering off a barrier, resonant behavior, proton transfer (quantum tunneling), and the classical/quantum correspondence principle. In subsequent papers, we will explore how the methods developed here may
be used to demonstrate such subjects. Computer-based instruction in quantum dynamics has been discussed in this Journal (1, 2). The value of using numerical solutions to the one-dimensional time-dependent Schrödinger equation in the physical chemistry course was discussed by Tanner (1). The existence of readily available and easily used software for numerically solving the timedependent Schrödinger equation, such as that discussed in this paper, should make it possible for instructors to use existing examples or to generate their own. We present here a specific implementation of a recently developed grid-based method that can be applied immediately by anyone having access to a spreadsheet. Using the template we describe, examples of the time-evolution of wavepackets in any potential can be generated easily. Given any initial wave function, the spreadsheet calculation determines the wave function propagated forward by the time step specified. Further time-evolution is carried out by using the propagated wave function as the initial wave function in a new calculation. When plotted, repeated propagation results in a sequence of graphs that can be interpreted as “frames” of a movie. We have chosen to implement these calculations on a spreadsheet for several reasons. First, spreadsheets are widely used and the methods described here should be portable to any common spreadsheet program. Thus, the existence of an easily used template should make wavepacket propagation accessible to everyone who wishes to use it. Spreadsheet templates have found similar utility in finding numerical solutions to the time-independent Schrödinger equation (3). Second, modern spreadsheet programs have sophisticated graphics tools making it relatively easy to prepare striking visual demonstrations of wavepacket motion. Finally, the recently developed computational techniques for wavepacket propagation described here lend themselves readily to this type of calculation. The next section of this paper discusses numerical wavepacket propagation in a coordinate-based grid representation, introducing the concept of the propagator matrix. It concludes by presenting an expression for the elements of a particular propagator matrix based upon distributed approximating functionals (DAFs). The implementation of this technique in the spreadsheet template is then described. Finally, the example of a wavepacket in a harmonic potential is demonstrated, with emphasis on how parameters such as the grid spacing and the time step are selected.
Vol. 74 No. 3 March 1997 • Journal of Chemical Education
335
Information • Textbooks • Media • Resources and
Wavepacket Propagation The propagation of a quantum mechanical wave function in time is governed by the time-dependent Schrödinger equation:
ih ∂ ψ x,t = Hψ x,t ∂t
(1)
ˆ is the Hamiltonian operator. For a very short where H time step τ we can always write that
ψ x,t + τ ≈ ψ x,t + τ
∂ψ x,t ∂t
(2)
(10)
V=V x
respectively. Equation 10 assumes that the potential energy depends only on the position of the particle. This is not always true (as, for example, in the case of a charged particle moving in an external magnetic field), but it is the common situation for a great many problems of chemical interest. In employing eq 4 as a computational tool we have to ˆ n to the wave function. This must apply operators such as H ˆ 2, be done with some care as can be seen by considering H that is
or from eq 1,
2
H = HH = K + V K + V ψ x,t + τ ≈ ψ x,t – iτ Hψ x,t h
(3)
In fact, this equation could be used for propagation if the time step were taken to be short enough. This equation can be made exact by adding contributions of higher powers of τ as follows:
ψ x,t + τ =
∞
Σ 1 n = 0 n!
ψ x,t + τ = P τ ψ x,t
ˆ is called the quantum mechaniwhere the operator P(τ) cal propagator and is given by ∞
Σ 1 n = 0 n!
–iHτ h
n
(6)
The sum in eq 6 is seen to be of the form of the Taylor series expansion of the exponential function, and indeed the propagator is commonly written in the form
P τ =exp –iHτ / h
(7)
In the template described here, the wave function is expressed in a so-called grid representation; that is, by its values on a set of evenly-spaced discrete points along the xaxis. If the wave function values on the grid are regarded to be a column matrix (vector), then a linear operator, like the propagator in eq 7, is represented by a square matrix. The problem becomes one of finding such a matrix representation of the propagator in eq 7 in a grid representation. ˆ can be written as The Hamiltonian H
H=K +V
(8)
where the kinetic and potential operators in the coordinate representation are given by 2 2 K =– h ∂ 2 2m ∂x
336
exp –iHτ / h ≠ exp –iKτ / h exp –iVτ / h
(4)
(5)
(9)
2
(11)
Note that H cannot be written in the usual binomial exˆ2 + 2K ˆ V ˆ +V ˆ 2) because K ˆ and V ˆ do pansion form (i.e., K not commute. As a result, it is easy to show that
n
–iHτ ψ x,t h
ˆ n times. That this is a ˆ n means application of H where H solution can be easily demonstrated by taking the derivative with respect to τ. It can be shown mathematically that this series solution converges for any value of τ. This means that no matter how large τ is, we can use eq 4 to calculate ψ(x,t + τ) from ψ(x,t), if we are willing to add enough terms. Thus eq 4 provides a computational means of calculating ψ(x,t + τ) if ψ(x,t ) is known. Eq 4 is often written in the operator form as
Pτ =
2
= K + KV + VK + V
(12)
as would be true for an ordinary exponential function. However, expressions of this kind are useful for short times even though they don’t hold exactly. A common approximation for the propagator of this kind is the socalled symmetric split-operator approximation (4):
exp –iHτ / h ≈ exp –iVτ / 2h exp –iKτ / h exp –iVτ / 2h
(13)
which is the most accurate split-operator approximation that involves only one application of the kinetic energy propagator. (Similar expressions can be derived from exact formal solutions to the time-dependent Schrödinger equation in which the splitting of the kinetic and potential energy is achieved using integrating factors (5, 6). These expressions can be used in place of the symmetric split operator in the spreadsheet without any additional computational effort.) If one expands the exponentials as Taylor series and collects terms, the τ0, τ1, and τ2 terms are exact, but the τn terms for n > 2 are only approximate. In this approximation, propagation is carried out by the successive application of three operators, the first and third being the potential depenˆ τ/2 h ). The second is exp({i K ˆ τ/ h ), dent operator, exp({iV the free propagator, so called because it represents the entire propagator in the absence of a potential. ˆ is Since according to eq 10, applying the operator V equivalent to multiplying by the potential, the result of apˆ τ/2 h ) to a function, ƒ(x), is simplying the operator exp({iV ˆ τ/2 h )f(x). In a grid representaply the new function exp({iV tion, in which a function is represented by a column matrix ˆ τ/2 h ) is repreof its values at the grid points {ƒ(xi)}, exp({iV sented as a simple diagonal matrix with its ith element equal to exp({iV(xi)τ/2 h ). ˆ (τ) = The application of the free propagator, F ˆ τ/ h ), to a function is not so simple. In general, the exp({i K result of applying an operator to a function must be written as an integral. Thus, to apply to a function we must evaluate the integral
Journal of Chemical Education • Vol. 74 No. 3 March 1997
∫ dx9F(x|x9;τ)ƒ(x9)
Information • Textbooks • Media • Resources not on those values further away. This implies that, owing to the finite amount of energy in the system, the particle can move, at most, only a certain distance in the time interval τ. The probability of propagating to a grid point i that is more than (w + 1) grid points away from the particular grid point xj in time τ is essentially zero. The bandedness greatly reduces the computational time required for matrix multiplication. Furthermore, the symmetry properties of this matrix require that only (w + 1) nonzero elements need to be calculated and stored. The spreadsheet described in this paper does wavepacket propagation using the split-operator formulation of eq 13, along with a representation of the free particle propagator that has the features described above. This recently developed propagator is based upon distributed approximating functionals (DAFs) (9–13). A detailed discussion of DAFs and the DAF propagator that involves some knowledge of Fourier transforms can be found in the references cited. This discussion is not absolutely essential for the use of the propagator spreadsheet, however. The important result is the DAF free propagator matrix, the discretized version of the propagator, which is given by eq 18, where σ is a width parameter that becomes complex for τ ≠ 0 with τ-dependence given by eq 19. The parameter ∆x is the spacing between grid points and H 2n is the Hermite polynomial of degree 2n. These polynomials can be calculated most efficiently using the standard recursion relations for Hermite polynomials. The DAF approximation involves truncating the sum as shown at a finite value of M, so that even Hermite polynomials up to order 2M are included. The spacing between grid points is ∆x. Because |xi – x j| = ∆x|i – j|, this propagator matrix displays the behavior discussed in the last section and we may express the matrix elements of F as in eq 17 with eq 20 in terms of p = |i – j|. The presence of the Gaussian factor exp[{(p∆x)2 / 2σ2(τ)] in eq 20 assures that ƒp(τ) is a rapidly decaying function of p and that may be assumed to be banded, since F ij = ƒ|i – j|. The expressions in eq 19 or 20 represent an approximation to the free propagator in terms of a Gaussian function multiplied by a sum of Hermite polynomials, both of complex argument. This approximation becomes exact in the limit as ∆x → 0 and M → ∞, regardless of the value of
where F(x|x9;τ) is called the matrix element of the free propagator. In a grid representation, the integral is approximated by numerical quadrature and correspondingly is carried out by discrete matrix multiplication:
dx′F x|x′;τ ƒ x′ ≈ Σ Fij τ ƒ x j
(14)
j
where F(τ) is the matrix representation of the free propagator on a coordinate grid. Using these results, eq 5 in the split operator approximation can be written as in eq 15, or on a discrete grid (eq 16). The evaluation of the integral in eq 15 or the matrix multiplication by F in eq 16 is the time-consuming part of the calculation. One way of doing the calculation is by applying a Fourier transform to obtain the wave function in the momentum representation in which the free-particle propagator is diagonal (4, 7, 8). After multiplying the diagonal free particle propagator and the momentum-space wavepacket, one Fourier transforms from momentum back to coordinate space. On a grid, this is a discrete Fourier transform and can be carried out using so-called “Fast Fourier Transforms”. Although the Fourier transform method has been used as a method for wavepacket propagation, better computational efficiency can be obtained if the free particle propagator is expressed directly in the coordinate representation. Under the appropriate conditions, it can be represented by a matrix having special properties which readily lend themselves to spreadsheet calculation. Specifically, the elements in this matrix depend only upon their distance from the diagonal (eq 17). That is, all the diagonal elements are equal to ƒ0 , all the elements adjacent to the diagonal are equal to ƒ1 , and so on. Furthermore, ƒp(τ) is a rapidly decreasing function of p = |i – j|, such that for greater than some number, w, the matrix elements are small enough to be assumed to be zero. Thus, the only nonzero elements of this matrix are concentrated in a band of width (2w + 1) around the diagonal. Such a matrix is said to be banded, with bandwidth (2w + 1). In physical terms, this means that the wave function propagated through a time step τ is dependent only on the previous values in a nearby region and
ψ x,t + τ = exp –iV x τ/2h
ψ x i,t + τ
=
dx′F x|x′;τ exp –iV x′ τ/2h ψ x′t
exp –iV x i τ/2h
Σj Fij τ
exp –iV x j τ/2h ψ x j t
Fij τ = ƒ|i– j| τ
Fij τ = ∆x
exp –(x i – x j ) 2 / 2σ 2(τ) 2π σ τ
M
Σ n=0
–1 4
exp – ( p∆x) 2 / 2σ 2(τ) 2π σ(τ)
M
Σ
n=0
–1 4
(16)
(17)
n
1 σ(0) n! σ(τ)
2n
H 2n
(x i – x j ) 2π σ τ
σ 2(τ) = σ 2(0) + ihτ / m
ƒ p = ∆x
(15)
n
1 σ(0) n! σ(τ)
(18)
(19)
2n
H 2n
p∆x 2 σ(τ)
(20)
Vol. 74 No. 3 March 1997 • Journal of Chemical Education
337
Information • Textbooks • Media • Resources σ(0), the initial (real) width parameter of the Gaussian. However, for finite M and ∆x > 0, there is an optimum value of σ(0), which depends upon both M and ∆x (11). With the stipulation that the optimum σ(0) is always used, the approximation is better for higher M. A general assessment of the errors inherent in this propagation scheme is somewhat complicated because of the various approximations that are made. These approximations include the use of the split-operator form of eq 13, the representation of the continuous variable x by a discrete grid, the truncation of the sum in eqs 19 and 20 at a finite value of M, and the restriction of the free propagator matrix to a finite bandwidth, (2w + 1). However, it is possible to calculate a few quantities from the propagator matrix elements and from the wave functions that allow one judge the goodness of the approximations used and to keep the errors within reasonable bounds. This is discussed more fully in the next section. The Propagation Spreadsheet Template This section describes in detail how the DAF propagation scheme on a grid, as defined by eqs 16–20, is implemented in a spreadsheet program. The program used was Microsoft Excel®, version 4.0 for the Macintosh. The entire calculation is carried out on one worksheet, which may be used as a template to generate other worksheets for specific examples. In addition to the worksheet, there are several macro programs which are contained in a separate macro sheet. The macro programs, while not strictly necessary, automate some of the most common anticipated manipulations of the worksheet. All of the information that is likely to be of interest to most users is contained in the first 8 columns of the template. Here, this will be referred to as the “public” section of the worksheet. Intermediate results and detailed calculations of the free propagator are in the so-called “private” section of the worksheet in the 9th column and beyond. In reality, there is no distinction between the two and the private section is readily accessed and altered by any user interested in the details of DAF propagation. The private section is simply not as conveniently placed as the public section. The primary data contained in the public section of the worksheet are the result of a single propagation calculation. This information is in a section of the worksheet labeled “Wave function table”. Each row of this table corresponds to one grid point and the table lists values of x and V(x) at that point as well as the values of the initial, or so-called “current” wave function and the wave function after it has been propagated forward in time by the specified time step. Because the wave functions are, in general, complex, two values are required for each of the wave function values, for the real and imaginary parts. The table also includes values of the squared magnitude (or probability function) for both the current and propagated wave functions. In addition to the wave function table, the public section also includes extensive directions explaining the use of the template for wavepacket propagation, which are intended to make it completely self-explanatory. There are also clearly marked cells for the user to enter values for various parameters necessary to define the particular propagation problem. Also included here are “buttons” associated with the programs on the macro sheet. Each macro can be executed by pressing the appropriate button. The worksheet calculates the propagated wave function in the wave function table from the current wave function using the split-operator formalism as shown in eq 16, with
338
a DAF representation for the free propagator matrix, as defined by eqs 17, 19, and 20. The parameter M is set to 10, thus including Hermite polynomials up to order 20 in the sum in eq 21. This value was chosen because it gave adequate results for all of the trial problems on which the template was used and did not require excessive computational time compared to the total time of carrying out a propagation step. Interested users should, however, find it relatively simple to change the value of M, as discussed below. The width parameter, σ(0), is set to the optimum value for M = 10. The elements of one band of the free propagator matrix are calculated from p = 0 to p = 20, thereby assuming a bandwidth less than or equal to (2 ? 20 + 1) = 41. The detailed implementation of this method within the worksheet is now described. It should be emphasized that, with the exception of the parameters defined by the user and the results of the propagation computation, the entire calculation takes place within the private section of the worksheet and need not concern the casual user. Six parameters dependent upon the specific calculation are contained in the appropriate marked cells in the public section of the worksheet. The parameters are the size of the time step for the propagator, called τ here but designated dt on the worksheet; the value of h ; the mass of the particle; the grid spacing, called ∆x here but designated dx on the worksheet; the minimum x on the grid; and the number of grid points. In the private section of the worksheet, the value of σ(0) is calculated in one cell from the formula σ(0) = 1.542∆x, which is the optimum value for M = 10. The value of σ(τ) is calculated from σ(0) using eq 19. Because this quantity is complex, real and imaginary parts must be separately calculated. For the convenience of subsequent computations based upon σ(τ), its magnitude and arguments as well as several other related quantities are also calculated in nearby cells. Cell references to σ(τ) are used in a table that computes values of the complex Gaussian function exp[ {(p∆x)2 / 2σ2 (τ)], for p = 0 to p = 20. These values are, in turn, used in a calculation of values of the complex Hermite functions, defined as the Gaussian function multiplied by the Hermite polynomial, or
HF(n, p) = exp –
( p∆x) 2 H 2n 2σ 2(τ)
p∆x 2 σ(τ)
(21)
These functions are required for the 21 values of p from 0 to 20, and the 11 values of n from 0 to 10. They are calculated in two tables, one for the real and another for the imaginary parts. The rows of each table correspond to values of n and the columns to values of p. The formulas in these tables make use of recursion relations for the Hermite functions so that each row beyond n = 2 uses identical formulas based upon the previous two rows. For this reason, it is a straightforward matter to extend the tables to higher values of n simply by copying the last row. The coefficients of the Hermite functions
C(n) =
∆x –1 2π σ(τ) 4
n
1 σ(0) n! σ(τ)
2n
(22)
are calculated in another table of complex numbers labeled by values of n. This table is similarly easy to extend to higher n values. Values of ƒp(τ) at each p are determined from the dot product of the coefficients and the values of the Hermite functions, thus:
Journal of Chemical Education • Vol. 74 No. 3 March 1997
Information • Textbooks • Media • Resources
ƒ p(τ) =
M
Σ C(n)HF(n, p) n=0
(23)
The formulas generating these values make use of an Excel function called “SUMPRODUCT,” which implements the dot product of two vectors. Finally, from the values of ƒp(τ) so calculated, a table of the 41 values in one row of the free propagator matrix is generated: {Fi{20,i;Fi{19,i;…;Fi,i;Fi+1,i ;…;F i+20,i } = {ƒ 20;ƒ 19;…;ƒ20} (24) making use of the symmetric nature of the free propagator. These matrix elements, which are independent of i, are contained in a table with 2 columns, for their real and imaginary parts, respectively. The validity of the free propagator so generated depends upon how banded the matrix is. Specifically, the bandwidth of the propagator must be less than or equal to 41 for the calculation to be valid because values outside of this band are not calculated. It turns out that the primary determining factor in the bandwidth of the free propagator matrix is the size of the time step. To assist the user in judging the validity of the approximations and deciding the value of τ, three quantities are generated from the free propagator. The first two quantities, labeled A and B on the template, are the sums of the real and imaginary parts, respectively, of the matrix elements of one row of F, from the table generated as described above. For an exact propagator, infinite in extent, these values should be 1 and 0, respectively. This can be seen most easily by considering the effect of the propagator upon a constant function. The timedependent Schrödinger equation shows that a constant function should be unchanged in time. This leads immediately to the requirement that the sum of one row of the free propagator matrix must be 1. Another quantity, labeled C on the template, is the sum of squares of the last 5 elements of the propagator row, from p = 16 to p = 20, divided by the sum of the squares of the 21 elements from p = 0 to p = 20: 20
C=
Σ p = 16
Fi + p,i
(25)
20
Σ
p=0
2
Fi + p,i
2
If the propagator matrix is sufficiently banded, the values of these matrix elements will have become very small for the elements at the end of the row and the value of C, which is positive definite, should be very small: C
