- Email: [email protected]
Solving the time-dependent Schroedinger equation by discarding high-energy basis functions
Chemical Physics Letters 501 (2010) 130–133
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Solving the time-dependent Schroedinger equation by discarding high-energy basis functions Mike McLeod, Tucker Carrington Jr. ⇑ Chemistry Department, Queen’s University, Kingston, Ontario, K7L 3N6 Canada
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 16 August 2010 In final form 15 October 2010 Available online 21 October 2010
It is demonstrated that it is possible to solve the time-dependent Schroedinger equation in 10d by discarding unimportant product basis functions. The matrix–vector products required to propagate an initial wave packet are done at a cost that scales as nDþ1 where n is a representative number of 1D basis functions and D is the dimensionality. Discarding unimportant basis functions drastically reduces the size of the vectors that must be stored. With D ¼ 20 and n ¼ 16 only 26 GB is required to store the wavefunction. The ideas are applied to a model Hamiltonian for hydrocarbon chains. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction To study the dynamics of a molecule or a dissociating system for which quantum effects are important one often solves the timedependent Schroedinger equation (TDSE) by representing the Hamiltonian and the evolving wave packet in a basis [1]. Numerical methods for solving the resulting matrix problem are well established [2,3]. The most versatile and effective methods require only the evaluation of matrix–vector products. Both the CPU and memory cost of such a calculation increase as the basis size increases. In this Letter we demonstrate that it is possible to solve a 10D TDSE to study energy flow in a hydrocarbon chain. Intramolecular vibrational energy transfer plays an important role in chemical reactions. It is frequently studied using classical molecular dynamics but quantum effects, for example zero-point energy, can be important. A natural way to evaluate a matrix–vector product is to make products of rows of the matrix with the vector. The memory cost of this approach scales as N 2 , where N is the size of the matrix. The most obvious basis is a direct product basis composed of nD products of functions of a single coordinate, where n is a representative number of basis functions for a single coordinate and D is the number of coordinates [4]. The increase in the size of a direct product basis with D is related to the so called curse of dimensionality [5]. For many interesting problems one needs about n ¼ 15 basis functions for each coordinate and therefore N 15D . A memory cost of 152D is prohibitive if D is larger than about 4. In addition, if D is large, such calculations require a lot of CPU time.
(M.
2. Model Hamiltonian The ideas are tested on a Hamiltonian that represents a linear hydrocarbon chain with one hydrogen atom at the end of a string of carbon atoms. It is similar to other Hamiltonians in the literature [12–15]. The Hamiltonian is,
bA þ H b C; b ¼H bH þ H H
⇑ Corresponding author. Fax: +1 613 533 6669. E-mail addresses: [email protected] queensu.ca (T. Carrington Jr.).
It is well-known that when a direct product basis is used it is possible to compute matrix–vector products at a cost that scales as nDþ1 , by doing summations sequentially [4,6–10]. This significantly reduces the CPU cost. Nonetheless, a direct product basis calculation with D = 10 is impossible because 4600 GB of memory are required to store a single (double precision) vector with 1510 components. In this Letter we show that calculations with D P 10 are possible. This is achieved by discarding unimportant basis functions without jeopardizing the favourable scaling of the CPU cost of the matrix–vector product. One way to reduce the computational cost of solving the Schroedinger equation is to use as basis functions products of specially adapted 1D functions that include the effect of coupling in an averaged sense. This is implemented in the multiconfiguration time-dependent Hartree approach [11]. The efficiency of this effective approach is related to its use of time-dependent basis functions, however, the multidimensional basis includes all the products of the 1D functions. In this Letter we propose an alternative strategy based on using only some multidimensional product functions.
McLeod),
Tucker.Carrington@
0009-2614/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2010.10.034
where
ð1Þ
M. McLeod, T. Carrington Jr. / Chemical Physics Letters 501 (2010) 130–133
bH ¼ H
L1 X p2i 1 þ lCC X2CC r2i ; 2 2 l CC i¼1
bA ¼ H
p2L
þ De ½1 earL 2 ; 2lCH L1 X pi piþ1 bC ¼ H : mc i¼1
ð2Þ ð3Þ ð4Þ
The C–C bond potentials are harmonic and the C–H bond potential is a Morse oscillator. The C–H bond is coupled to the adjacent C–C by a kinetic coupling term; mc is the mass of the carbon atom. The dissociation energy of the Morse oscillator is De; lCH and lCC are reduced masses, XCC is the C–C frequency, and De , a are the usual Morse parameters. L is the length of the chain. The harmonic frequency is XCC = 0.008429 a.u., the Morse parameters are De = 0.134 a.u. and a ¼ 1:03a1 0 . Potential and other coupling terms could be included. We want to choose 1D basis functions so that some of the product basis functions 1 2 i1 ðr 1 Þ i2 ðr 2 Þ
Ui1 ;i2 ;...;iD ¼ a
a
D iD ðr L Þ;
a
ð5Þ
can be discarded (see the next section). One way to do this is to use akik ðrk Þ functions that approximate solutions of 1D Hamiltonians extracted from the full Hamiltonian by neglecting coupling. It is natural to use harmonic oscillator basis functions for the C–C bonds. For the C–H bond we could use Morse basis functions. However, doing so would complicate the computation of matrix elements of b C . Instead, we use a transformation suggested the i ¼ L 1 term in H by Cooper [16] and use harmonic basis functions. First, we re-write the Hamiltonian in terms of dimensionless coordinates
b b wn þ H b wn ¼ H ¼ H b wn þ H b wn þ H b wn ; H H A C1 C2 hc
ð6Þ
In terms of Q L and P L
1 1 wn d H xCH P^n 2 þ Q 2n bxm ½P^n ; ½P^n ; Q n þ þ A ¼ 2 4 1 2 2 c þ b xm ½ P n ; Q n þ ; 8
ð20Þ
where ½Pn ; Q n þ is an anticommutator. Written in creation and annib wn becomes, hilation operators H A wn d d d d H A ¼ HA1 þ H A2 þ HA3 ;
ð21Þ
where
cþ c 1 ; d H A1 ¼ xm hc Aa Aa þ 2
ð22Þ
n o 1 3 d cþ c cþ c cþ d c A d H Aþa 3 A ; A2 ¼ pffiffiffi bxm hc a a Aa Aa Aa Aa þ Aa 2 2
ð23Þ
1 2 d H A3 ¼ b xm 8
n þ 2 d 2 d 4o d Aa 2 d Aa Aa Aa : Aþa 4 d Aþa 2 þ d ð24Þ
To compute matrix elements we write the final Hamiltonian in terms of raising and lowering operators. 3. Propagation method We propagate using the Chebyshev method [17,18]. To determine the evolving wave packet at intermediate times one can divide the total propagation interval into steps and write
eiHDt ¼ eiðcMþdIÞ ;
ð25Þ
where H and M are matrices. In this equation,
where L1 X xCC 2 Pi þ Q 2i ; 2 i¼1 xCH ~2 ¼ pL þ Q 2L ; 2 L2 X lCC hXCC Pi Piþ1 ¼ ; mc i¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~L PL1 ahp ¼ lCC hXCC ; bmc
131
c ¼ ðDt=2Þðkmax kmin Þ;
d ¼ ðDt=2Þðkmax þ kmin Þ;
ð26Þ
b wn ¼ H H
ð7Þ
and
b wn H A
ð8Þ
M¼
ð9Þ
with kmax and kmin being the largest and smallest eigenvalues of H. A Chebyshev series is used to expand eicM ,
b wn H C1 b wn H C2
ð10Þ
H Iðkmax þ kmin Þ ; kmax kmin
eicM ¼
ð27Þ
1 X ð2 dk;0 ÞðiÞk J k ðcðDtÞÞT k ðMÞ:
ð28Þ
k¼0
with
1 ebqL QL ¼ ; b sffiffiffiffiffiffiffiffiffi 2x ; b¼
ð11Þ
AðtÞ ¼ hwð0Þ j wðtÞi ¼ eid ð12Þ
xCH
ð13Þ ð14Þ ð15Þ ð16Þ
lCH
1 Pi ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pi ; lCC hXCC rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
xCC
lCC XCC
h XCC h ¼ ; hc
ri ;
1 X ð2 dk;0 ÞðiÞk J k ðcðDtÞÞhwð0Þ j wk i; k¼0
ð29Þ
a qL ¼ r L ; b b ~L ¼ p; p ah L 2 a2 h ; x¼ 2hclCH sffiffiffiffiffiffiffiffiffiffiffiffiffi 2a2 De ; xCH hc ¼ h
Qi ¼
To compute a survival probability there is no need to take steps. Instead one can use [18]
i ¼ 1; . . . ; L 1; i ¼ 1; . . . ; L 1;
i ¼ 1; . . . ; L 1:
ð17Þ ð18Þ ð19Þ
where j wk i ¼ T k ðMÞ j wð0Þi and compute a survival probability PðtÞ from j AðtÞj2 . Rather than using the full direct product basis we use part of it. The resulting basis is a nondirect product basis, but selected from a direct product basis. The idea of using only some of the functions of a direct product basis has been previously exploited [19–26]. The functions we retain are selected by constraining the indices in the summations done to compute matrix–vector products. We exclude basis functions associated with large zeroth-order energies. This can be done without sacrificing the favourable nDþ1 scaling relation for the cost of a single matrix–vector product. The same idea was used to compute energy levels in Ref. [26]. If zeroth-order energies for all the coordinates are similar, and if coupling is equally important for all the coordinates, one would exclude functions for which the sum of the 1D quantum numbers
132
M. McLeod, T. Carrington Jr. / Chemical Physics Letters 501 (2010) 130–133
is larger than some chosen value, i.e., exclude all functions Ui1 ;i2 ;...;iD for which i1 þ i2 þ þ iD > N t , where N t is a threshold value and ik ¼ 0; 1; 2; . . . ; N t . This sort of constraint greatly reduces the size of the basis. For example, for a Hamiltonian with D ¼ 10 and N t ¼ 15, the size of the direct product basis is 1610 1 1012 but t Þ! the size of the constrained basis is ðDþN 3 106 . When D ¼ 15 D!N t ! the size of the basis is reduced from 1:2 1018 to 1:6 108 . In the hydrocarbon chain problem the C–H coordinate does not have the same zeroth-order energy as the C–C coordinates and is more strongly coupled. We therefore exclude basis functions for which i1 þ i2 þ þ iL1 > N t and keep N t þ 1 values of iL . It is probably possible to further reduce the size of the basis but our reduction scheme is good enough to demonstrate the ideas. Excluding unnecessary functions from the basis: (1) reduces the size of the vectors one must store and therefore makes it possible to do calculations with D P 10; (2) accelerates the solution of the TDSE because it reduces the cost of evaluating a single matrix–vector product and also because it decreases the spectral range and hence the number of matrix–vector products required. Matrix–vector products are done term by term. The matrix representation of each term in the Hamiltonian can be written as a product of small matrices ðF1Þi1 ;i 0 ðF2Þi2 ;i 0 ðFDÞiD ;iD0 . At most two 1 2 of the factors are not identity matrices. For each of the terms in wn b only one factor is not an identity and it is diagonal. For a H H b wn , for which the i5 and i6 fachydrocarbon chain with L ¼ 6 and H C2 tors are not identity matrices, the matrix vector product is evaluated as
wi1 ;i2 ;i3 ;i4 ;i5 ;i6 ¼
N t i1X i2 i3 i4 i05
ðF5Þi5 ;i05
Nt X ðF6Þi6 ;i0 v i1 ;i2 ;i3 ;i4 ;i05 ;i06 ; 6
i06
ð30Þ
where underlined indices are subject to the constraint that their b wn term for which the i4 sum is less than or equal to N t . For a H C1 and i5 factors are not identity matrices the matrix vector product is evaluated as
wi1 ;i2 ;i3 ;i4 ;i5 ;i6 ¼
N t iX 1 i2 i3 i04
ðF4Þi4 ;i0
4
Nt i1X i2 i3 i4 i05
ðF5Þi5 ;i05 v i1 ;i2 ;i3 ;i04 ;i05 ;i6 :
Figure 2. L ¼ 10; N t ¼ 13 (circles) and 14 (squares); Reduced basis sizes: 6 963 880 and 12 257 850 (direct product bases would have 6 1011 and 1 1012 functions).
tion corresponding to the v ¼ 5 state of the C–H Morse oscillator and 1D ground state wavefunctions for each of the C–C coordinates. The Morse v ¼ 5 wavefunction is determined by solving b A in the harmonic basis used for the for the eigenfunctions of H C–H coordinate. In Figure 1 we show the survival probability for L ¼ 6. We have verified that we obtain very nearly the same results with the constrained and full direct product bases (see Figure 1). This calculation could be done with a direct product basis. Our results agree well with those of Ref. [15] for the same chains. To clearly demonstrate the utility of constraining the basis. We have also done calculations for L ¼ 10 (the maximum chain length in Ref. [15] is 6). The results are shown in Figure 2. In the L ¼ 10 case propagating to 35 fs takes 11 h on a modern desktop personal computer.
ð31Þ
4. Results For L = 2, 6, and 10, (i.e. 2, 6, 10 carbon atoms and one hydrogen atom) we have computed the survival probability, PðtÞ ¼ jhwð0ÞjwðtÞij2 , where wð0Þ is a product of the 1D wavefunc-
Figure 1. L ¼ 6; N t ¼ 13 (squares) and 14 (circles). Reduced basis sizes: 119 952 and 174 420. Direct product basis results with N t ¼ 13 are triangles.
5. Conclusion In this Letter we have demonstrated that by using only part of a direct product basis it is possible to solve the TDSE to study intramolecular vibrational energy transfer in 10D. The model Hamiltonian we use is simple but realistic and has been used with other computational approaches to study energy flow in hydrocarbon chains. Our basis reduction scheme allows us to obtain accurate results. It has two key advantages. First, the basis size is small. For a problem with similar zeroth-order 1D energies and n ¼ 15 (number of basis functions for a single coordinate) the size of a direct product basis is: 5 1011 , for D ¼ 10; 4 1017 , for D ¼ 15; 3 1023 , for D ¼ 20. Even with D ¼ 10 it would take 4000 GB to store one vector. On the other hand, if basis functions with i1 þ i2 þ þ iD > N t ¼ n 1 can be discarded and N t ¼ 14 the size of the product basis is: 2:0 106 , for D ¼ 10; 7:7 107 , for D ¼ 15; 1:4 109 , for D ¼ 20. Even the D ¼ 20 calculation requires only 11 GB! The basis size increase is less than linear. The Hamiltonian we use has only bilinear momentum coupling terms. Potential terms that couple neighbouring coordinates might also be important but terms coupling coordinates that are well separated are probably unimportant. Therefore, simple Hamiltonians with a fairly small number of terms should be realistic enough for many energy transfer problems. Potential (and additional kinetic) coupling terms could be treated in the same way we have treated the bilinear terms. If the number of potential terms were large it would be better to use quadrature. To avoid using a quadrature grid much larger than the basis one would need to
M. McLeod, T. Carrington Jr. / Chemical Physics Letters 501 (2010) 130–133
use Smolyak-based ideas similar to those in Ref. [27]. By treating part of the problem exactly, explicitly, and quantum mechanically and the rest approximately, implicitly or classically, ideas presented in this Letter could be used with a system-bath approach for even larger problems. Acknowledgements Calculations were done on a computer bought with funds from the Canadian Foundation for Innovation. This work was supported by the Natural Sciences and Engineering Research Council of Canada. References [1] [2] [3] [4] [5]
John Z.H. Zhang, World Scientific Publishing Company (1999). R. Kosloff, Annu. Rev. Chem. Phys. 45 (1994) 145. C. Leforestier et al., J. Comp. Phys. 94 (1991) 59. J.C. Light, T. Carrington Jr., Adv. Chem. Phys. 114 (2000) 263. D. Donoho, High-dimensional data analysis: the curses and blessings of dimensionality, mathematical challenges of the 21st century (2000).
133
[6] M.J. Bramley, T. Carrington, J. Chem. Phys. 99 (1993) 8519. [7] X.G. Wang, T. Carrington, J. Chem. Phys. 118 (2002) 6946. [8] T. Carrington Jr., in: Paul von Ragué Schleyer (Ed.), Encyclopedia of Computational Chemistry, vol. 5, John Wiley & Sons, 1998. [9] X.G. Wang, T. Carrington, J. Chem. Phys. 115 (2001) 9781. [10] Rongqing Chen, Guobin Ma, Hua Guo, J. Chem. Phys. 114 (2001) 4763. [11] M.H. Beck, A. Jaeckle, G.A. Worth, H.-D. Meyer, Phys. Rep. 324 (2000) 1. [12] Alessandro Lami, Giovanni Villani, J. Chem. Phys. 90 (1989) 3559. [13] S. Mitra, S.D. Schwartz, J. Chem. Phys. 104 (1995) 7539. [14] D. Gelman, S.D. Schwartz, J. Chem. Phys. 130 (2009) 134110. [15] Maria Topaler, Nancy Makri, J. Chem. Phys. 97 (1992) 9001. [16] I.L. Cooper, Chem. Phys. 112 (1987) 67. [17] H. Tal-Ezer, R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [18] R. Kosloff, J. Phys. Chem. 92 (1988) 2087. [19] R. Dawes, T. Carrington, J. Chem. Phys. 122 (2005) 134101. [20] L. Halonen, M.S. Child, J. Chem. Phys. 79 (1983) 4355. [21] J.M. Bowman, S. Carter, X. Huang, Int. Rev. Phys. Chem. 22 (2003) 533. [22] Didier Bégué, Neil Gohaud, Claude Pouchan, Patrick Cassam-Chenai, Jacques Liévin, J. Chem. Phys. 127 (2007) 164115. [23] Andrew Maynard, Robert E. Wyatt, Christopher Iung, J. Chem. Phys. 106 (1997) 9483. [24] B. Poirier, J. Theo. Comput. Chem. 2 (2003) 65. [25] Jason Cooper, Tucker Carrington, J. Chem. Phys. 130 (2009) 214110. [26] X.-G. Wang, T. Carrington Jr., J. Phys. Chem. A 105 (2001) 2575. [27] Gustavo Avila, Tucker Carrington, J. Chem. Phys. 131 (2009) 174103-1.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Solving the time-dependent Schroedinger equation by discarding high-energy basis functions Mike McLeod, Tucker Carrington Jr. ⇑ Chemistry Department, Queen’s University, Kingston, Ontario, K7L 3N6 Canada
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 16 August 2010 In final form 15 October 2010 Available online 21 October 2010
It is demonstrated that it is possible to solve the time-dependent Schroedinger equation in 10d by discarding unimportant product basis functions. The matrix–vector products required to propagate an initial wave packet are done at a cost that scales as nDþ1 where n is a representative number of 1D basis functions and D is the dimensionality. Discarding unimportant basis functions drastically reduces the size of the vectors that must be stored. With D ¼ 20 and n ¼ 16 only 26 GB is required to store the wavefunction. The ideas are applied to a model Hamiltonian for hydrocarbon chains. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction To study the dynamics of a molecule or a dissociating system for which quantum effects are important one often solves the timedependent Schroedinger equation (TDSE) by representing the Hamiltonian and the evolving wave packet in a basis [1]. Numerical methods for solving the resulting matrix problem are well established [2,3]. The most versatile and effective methods require only the evaluation of matrix–vector products. Both the CPU and memory cost of such a calculation increase as the basis size increases. In this Letter we demonstrate that it is possible to solve a 10D TDSE to study energy flow in a hydrocarbon chain. Intramolecular vibrational energy transfer plays an important role in chemical reactions. It is frequently studied using classical molecular dynamics but quantum effects, for example zero-point energy, can be important. A natural way to evaluate a matrix–vector product is to make products of rows of the matrix with the vector. The memory cost of this approach scales as N 2 , where N is the size of the matrix. The most obvious basis is a direct product basis composed of nD products of functions of a single coordinate, where n is a representative number of basis functions for a single coordinate and D is the number of coordinates [4]. The increase in the size of a direct product basis with D is related to the so called curse of dimensionality [5]. For many interesting problems one needs about n ¼ 15 basis functions for each coordinate and therefore N 15D . A memory cost of 152D is prohibitive if D is larger than about 4. In addition, if D is large, such calculations require a lot of CPU time.
(M.
2. Model Hamiltonian The ideas are tested on a Hamiltonian that represents a linear hydrocarbon chain with one hydrogen atom at the end of a string of carbon atoms. It is similar to other Hamiltonians in the literature [12–15]. The Hamiltonian is,
bA þ H b C; b ¼H bH þ H H
⇑ Corresponding author. Fax: +1 613 533 6669. E-mail addresses: [email protected] queensu.ca (T. Carrington Jr.).
It is well-known that when a direct product basis is used it is possible to compute matrix–vector products at a cost that scales as nDþ1 , by doing summations sequentially [4,6–10]. This significantly reduces the CPU cost. Nonetheless, a direct product basis calculation with D = 10 is impossible because 4600 GB of memory are required to store a single (double precision) vector with 1510 components. In this Letter we show that calculations with D P 10 are possible. This is achieved by discarding unimportant basis functions without jeopardizing the favourable scaling of the CPU cost of the matrix–vector product. One way to reduce the computational cost of solving the Schroedinger equation is to use as basis functions products of specially adapted 1D functions that include the effect of coupling in an averaged sense. This is implemented in the multiconfiguration time-dependent Hartree approach [11]. The efficiency of this effective approach is related to its use of time-dependent basis functions, however, the multidimensional basis includes all the products of the 1D functions. In this Letter we propose an alternative strategy based on using only some multidimensional product functions.
McLeod),
Tucker.Carrington@
0009-2614/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2010.10.034
where
ð1Þ
M. McLeod, T. Carrington Jr. / Chemical Physics Letters 501 (2010) 130–133
bH ¼ H
L1 X p2i 1 þ lCC X2CC r2i ; 2 2 l CC i¼1
bA ¼ H
p2L
þ De ½1 earL 2 ; 2lCH L1 X pi piþ1 bC ¼ H : mc i¼1
ð2Þ ð3Þ ð4Þ
The C–C bond potentials are harmonic and the C–H bond potential is a Morse oscillator. The C–H bond is coupled to the adjacent C–C by a kinetic coupling term; mc is the mass of the carbon atom. The dissociation energy of the Morse oscillator is De; lCH and lCC are reduced masses, XCC is the C–C frequency, and De , a are the usual Morse parameters. L is the length of the chain. The harmonic frequency is XCC = 0.008429 a.u., the Morse parameters are De = 0.134 a.u. and a ¼ 1:03a1 0 . Potential and other coupling terms could be included. We want to choose 1D basis functions so that some of the product basis functions 1 2 i1 ðr 1 Þ i2 ðr 2 Þ
Ui1 ;i2 ;...;iD ¼ a
a
D iD ðr L Þ;
a
ð5Þ
can be discarded (see the next section). One way to do this is to use akik ðrk Þ functions that approximate solutions of 1D Hamiltonians extracted from the full Hamiltonian by neglecting coupling. It is natural to use harmonic oscillator basis functions for the C–C bonds. For the C–H bond we could use Morse basis functions. However, doing so would complicate the computation of matrix elements of b C . Instead, we use a transformation suggested the i ¼ L 1 term in H by Cooper [16] and use harmonic basis functions. First, we re-write the Hamiltonian in terms of dimensionless coordinates
b b wn þ H b wn ¼ H ¼ H b wn þ H b wn þ H b wn ; H H A C1 C2 hc
ð6Þ
In terms of Q L and P L
1 1 wn d H xCH P^n 2 þ Q 2n bxm ½P^n ; ½P^n ; Q n þ þ A ¼ 2 4 1 2 2 c þ b xm ½ P n ; Q n þ ; 8
ð20Þ
where ½Pn ; Q n þ is an anticommutator. Written in creation and annib wn becomes, hilation operators H A wn d d d d H A ¼ HA1 þ H A2 þ HA3 ;
ð21Þ
where
cþ c 1 ; d H A1 ¼ xm hc Aa Aa þ 2
ð22Þ
n o 1 3 d cþ c cþ c cþ d c A d H Aþa 3 A ; A2 ¼ pffiffiffi bxm hc a a Aa Aa Aa Aa þ Aa 2 2
ð23Þ
1 2 d H A3 ¼ b xm 8
n þ 2 d 2 d 4o d Aa 2 d Aa Aa Aa : Aþa 4 d Aþa 2 þ d ð24Þ
To compute matrix elements we write the final Hamiltonian in terms of raising and lowering operators. 3. Propagation method We propagate using the Chebyshev method [17,18]. To determine the evolving wave packet at intermediate times one can divide the total propagation interval into steps and write
eiHDt ¼ eiðcMþdIÞ ;
ð25Þ
where H and M are matrices. In this equation,
where L1 X xCC 2 Pi þ Q 2i ; 2 i¼1 xCH ~2 ¼ pL þ Q 2L ; 2 L2 X lCC hXCC Pi Piþ1 ¼ ; mc i¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~L PL1 ahp ¼ lCC hXCC ; bmc
131
c ¼ ðDt=2Þðkmax kmin Þ;
d ¼ ðDt=2Þðkmax þ kmin Þ;
ð26Þ
b wn ¼ H H
ð7Þ
and
b wn H A
ð8Þ
M¼
ð9Þ
with kmax and kmin being the largest and smallest eigenvalues of H. A Chebyshev series is used to expand eicM ,
b wn H C1 b wn H C2
ð10Þ
H Iðkmax þ kmin Þ ; kmax kmin
eicM ¼
ð27Þ
1 X ð2 dk;0 ÞðiÞk J k ðcðDtÞÞT k ðMÞ:
ð28Þ
k¼0
with
1 ebqL QL ¼ ; b sffiffiffiffiffiffiffiffiffi 2x ; b¼
ð11Þ
AðtÞ ¼ hwð0Þ j wðtÞi ¼ eid ð12Þ
xCH
ð13Þ ð14Þ ð15Þ ð16Þ
lCH
1 Pi ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pi ; lCC hXCC rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
xCC
lCC XCC
h XCC h ¼ ; hc
ri ;
1 X ð2 dk;0 ÞðiÞk J k ðcðDtÞÞhwð0Þ j wk i; k¼0
ð29Þ
a qL ¼ r L ; b b ~L ¼ p; p ah L 2 a2 h ; x¼ 2hclCH sffiffiffiffiffiffiffiffiffiffiffiffiffi 2a2 De ; xCH hc ¼ h
Qi ¼
To compute a survival probability there is no need to take steps. Instead one can use [18]
i ¼ 1; . . . ; L 1; i ¼ 1; . . . ; L 1;
i ¼ 1; . . . ; L 1:
ð17Þ ð18Þ ð19Þ
where j wk i ¼ T k ðMÞ j wð0Þi and compute a survival probability PðtÞ from j AðtÞj2 . Rather than using the full direct product basis we use part of it. The resulting basis is a nondirect product basis, but selected from a direct product basis. The idea of using only some of the functions of a direct product basis has been previously exploited [19–26]. The functions we retain are selected by constraining the indices in the summations done to compute matrix–vector products. We exclude basis functions associated with large zeroth-order energies. This can be done without sacrificing the favourable nDþ1 scaling relation for the cost of a single matrix–vector product. The same idea was used to compute energy levels in Ref. [26]. If zeroth-order energies for all the coordinates are similar, and if coupling is equally important for all the coordinates, one would exclude functions for which the sum of the 1D quantum numbers
132
M. McLeod, T. Carrington Jr. / Chemical Physics Letters 501 (2010) 130–133
is larger than some chosen value, i.e., exclude all functions Ui1 ;i2 ;...;iD for which i1 þ i2 þ þ iD > N t , where N t is a threshold value and ik ¼ 0; 1; 2; . . . ; N t . This sort of constraint greatly reduces the size of the basis. For example, for a Hamiltonian with D ¼ 10 and N t ¼ 15, the size of the direct product basis is 1610 1 1012 but t Þ! the size of the constrained basis is ðDþN 3 106 . When D ¼ 15 D!N t ! the size of the basis is reduced from 1:2 1018 to 1:6 108 . In the hydrocarbon chain problem the C–H coordinate does not have the same zeroth-order energy as the C–C coordinates and is more strongly coupled. We therefore exclude basis functions for which i1 þ i2 þ þ iL1 > N t and keep N t þ 1 values of iL . It is probably possible to further reduce the size of the basis but our reduction scheme is good enough to demonstrate the ideas. Excluding unnecessary functions from the basis: (1) reduces the size of the vectors one must store and therefore makes it possible to do calculations with D P 10; (2) accelerates the solution of the TDSE because it reduces the cost of evaluating a single matrix–vector product and also because it decreases the spectral range and hence the number of matrix–vector products required. Matrix–vector products are done term by term. The matrix representation of each term in the Hamiltonian can be written as a product of small matrices ðF1Þi1 ;i 0 ðF2Þi2 ;i 0 ðFDÞiD ;iD0 . At most two 1 2 of the factors are not identity matrices. For each of the terms in wn b only one factor is not an identity and it is diagonal. For a H H b wn , for which the i5 and i6 fachydrocarbon chain with L ¼ 6 and H C2 tors are not identity matrices, the matrix vector product is evaluated as
wi1 ;i2 ;i3 ;i4 ;i5 ;i6 ¼
N t i1X i2 i3 i4 i05
ðF5Þi5 ;i05
Nt X ðF6Þi6 ;i0 v i1 ;i2 ;i3 ;i4 ;i05 ;i06 ; 6
i06
ð30Þ
where underlined indices are subject to the constraint that their b wn term for which the i4 sum is less than or equal to N t . For a H C1 and i5 factors are not identity matrices the matrix vector product is evaluated as
wi1 ;i2 ;i3 ;i4 ;i5 ;i6 ¼
N t iX 1 i2 i3 i04
ðF4Þi4 ;i0
4
Nt i1X i2 i3 i4 i05
ðF5Þi5 ;i05 v i1 ;i2 ;i3 ;i04 ;i05 ;i6 :
Figure 2. L ¼ 10; N t ¼ 13 (circles) and 14 (squares); Reduced basis sizes: 6 963 880 and 12 257 850 (direct product bases would have 6 1011 and 1 1012 functions).
tion corresponding to the v ¼ 5 state of the C–H Morse oscillator and 1D ground state wavefunctions for each of the C–C coordinates. The Morse v ¼ 5 wavefunction is determined by solving b A in the harmonic basis used for the for the eigenfunctions of H C–H coordinate. In Figure 1 we show the survival probability for L ¼ 6. We have verified that we obtain very nearly the same results with the constrained and full direct product bases (see Figure 1). This calculation could be done with a direct product basis. Our results agree well with those of Ref. [15] for the same chains. To clearly demonstrate the utility of constraining the basis. We have also done calculations for L ¼ 10 (the maximum chain length in Ref. [15] is 6). The results are shown in Figure 2. In the L ¼ 10 case propagating to 35 fs takes 11 h on a modern desktop personal computer.
ð31Þ
4. Results For L = 2, 6, and 10, (i.e. 2, 6, 10 carbon atoms and one hydrogen atom) we have computed the survival probability, PðtÞ ¼ jhwð0ÞjwðtÞij2 , where wð0Þ is a product of the 1D wavefunc-
Figure 1. L ¼ 6; N t ¼ 13 (squares) and 14 (circles). Reduced basis sizes: 119 952 and 174 420. Direct product basis results with N t ¼ 13 are triangles.
5. Conclusion In this Letter we have demonstrated that by using only part of a direct product basis it is possible to solve the TDSE to study intramolecular vibrational energy transfer in 10D. The model Hamiltonian we use is simple but realistic and has been used with other computational approaches to study energy flow in hydrocarbon chains. Our basis reduction scheme allows us to obtain accurate results. It has two key advantages. First, the basis size is small. For a problem with similar zeroth-order 1D energies and n ¼ 15 (number of basis functions for a single coordinate) the size of a direct product basis is: 5 1011 , for D ¼ 10; 4 1017 , for D ¼ 15; 3 1023 , for D ¼ 20. Even with D ¼ 10 it would take 4000 GB to store one vector. On the other hand, if basis functions with i1 þ i2 þ þ iD > N t ¼ n 1 can be discarded and N t ¼ 14 the size of the product basis is: 2:0 106 , for D ¼ 10; 7:7 107 , for D ¼ 15; 1:4 109 , for D ¼ 20. Even the D ¼ 20 calculation requires only 11 GB! The basis size increase is less than linear. The Hamiltonian we use has only bilinear momentum coupling terms. Potential terms that couple neighbouring coordinates might also be important but terms coupling coordinates that are well separated are probably unimportant. Therefore, simple Hamiltonians with a fairly small number of terms should be realistic enough for many energy transfer problems. Potential (and additional kinetic) coupling terms could be treated in the same way we have treated the bilinear terms. If the number of potential terms were large it would be better to use quadrature. To avoid using a quadrature grid much larger than the basis one would need to
M. McLeod, T. Carrington Jr. / Chemical Physics Letters 501 (2010) 130–133
use Smolyak-based ideas similar to those in Ref. [27]. By treating part of the problem exactly, explicitly, and quantum mechanically and the rest approximately, implicitly or classically, ideas presented in this Letter could be used with a system-bath approach for even larger problems. Acknowledgements Calculations were done on a computer bought with funds from the Canadian Foundation for Innovation. This work was supported by the Natural Sciences and Engineering Research Council of Canada. References [1] [2] [3] [4] [5]
John Z.H. Zhang, World Scientific Publishing Company (1999). R. Kosloff, Annu. Rev. Chem. Phys. 45 (1994) 145. C. Leforestier et al., J. Comp. Phys. 94 (1991) 59. J.C. Light, T. Carrington Jr., Adv. Chem. Phys. 114 (2000) 263. D. Donoho, High-dimensional data analysis: the curses and blessings of dimensionality, mathematical challenges of the 21st century (2000).
133
[6] M.J. Bramley, T. Carrington, J. Chem. Phys. 99 (1993) 8519. [7] X.G. Wang, T. Carrington, J. Chem. Phys. 118 (2002) 6946. [8] T. Carrington Jr., in: Paul von Ragué Schleyer (Ed.), Encyclopedia of Computational Chemistry, vol. 5, John Wiley & Sons, 1998. [9] X.G. Wang, T. Carrington, J. Chem. Phys. 115 (2001) 9781. [10] Rongqing Chen, Guobin Ma, Hua Guo, J. Chem. Phys. 114 (2001) 4763. [11] M.H. Beck, A. Jaeckle, G.A. Worth, H.-D. Meyer, Phys. Rep. 324 (2000) 1. [12] Alessandro Lami, Giovanni Villani, J. Chem. Phys. 90 (1989) 3559. [13] S. Mitra, S.D. Schwartz, J. Chem. Phys. 104 (1995) 7539. [14] D. Gelman, S.D. Schwartz, J. Chem. Phys. 130 (2009) 134110. [15] Maria Topaler, Nancy Makri, J. Chem. Phys. 97 (1992) 9001. [16] I.L. Cooper, Chem. Phys. 112 (1987) 67. [17] H. Tal-Ezer, R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [18] R. Kosloff, J. Phys. Chem. 92 (1988) 2087. [19] R. Dawes, T. Carrington, J. Chem. Phys. 122 (2005) 134101. [20] L. Halonen, M.S. Child, J. Chem. Phys. 79 (1983) 4355. [21] J.M. Bowman, S. Carter, X. Huang, Int. Rev. Phys. Chem. 22 (2003) 533. [22] Didier Bégué, Neil Gohaud, Claude Pouchan, Patrick Cassam-Chenai, Jacques Liévin, J. Chem. Phys. 127 (2007) 164115. [23] Andrew Maynard, Robert E. Wyatt, Christopher Iung, J. Chem. Phys. 106 (1997) 9483. [24] B. Poirier, J. Theo. Comput. Chem. 2 (2003) 65. [25] Jason Cooper, Tucker Carrington, J. Chem. Phys. 130 (2009) 214110. [26] X.-G. Wang, T. Carrington Jr., J. Phys. Chem. A 105 (2001) 2575. [27] Gustavo Avila, Tucker Carrington, J. Chem. Phys. 131 (2009) 174103-1.