Computation of interior Eigenstates of Large Matrices using the quasiadiabatic evolution of

instantaneous eigenvectors

Georges JOLICARD

Observatoire de Besançon. Université de Franche Comté, 41, Av. de l’Observatoire, BP 1615, 25010 Besançon Cedex, France

John P. KILLINGBECK

department of Applied Mathematics, The University, Cottingham Road HULL, HU6 7RX, United Kingdom

February 14, 1995

Abstract

A two-stage iterative scheme is proposed to handle a central problem of molecular dynamics, the computation of interior eigenvalues of large Hamiltonian matrices. The proposed method involves an initial propagation process for a time-dependent wave operator which is then inserted in an iterative process (R.D.W.A. or S.C.M.) to yield the exact stationary wave operator. The merits of the wave operator formalism for quasiadiabatic propagation are analysed, and possible improvements such as the use of partial adiabatic representations and spectral filters, are outlined. The proposed algorithm is applied to the test case of two coupled oscillators with variable coupling strength, and yields accurate results even with small switching times.

1 INTRODUCTION

Numerous theoretical studies in molecular and chemical physics use discrete representations of the molecular continua and, via the Floquet treatment of the field-matter interaction, of the time variable. These Discrete Variable Representations (D.V.R.) lead to complex eigenvalue problems. Well known examples include the study of the molecular resonance states which influence the vibrational energy redistribution and fragmentation of Van der Waals molecules and the calculation of the resonances induced by the field-matter interaction in photodissociation processes. [1, 2, 3]

In most such cases the eigenproblems which arise are character-

ized by several common features : The matrix representation of the Hamiltonian is very large (typically N > 10^{4}); these matrice are complex, as a consequence of absorbing boundary conditions introduced by similarity transformations in the complex plane [4, 5] or by the addition of an optical potential [6]; the number of eigenvectors effectively contributing to the fragmentation processes studied is often small and only a small set of eigenstates is needed.

Complete direct diagonalization of the matrices concerned is

out of the question and over the last decade a large variety of methods has been proposed to solve these eigenproblems. These methods fall into two groups: indirect and direct methods. The indirect methods involve wave packet propagation and extract the required results from some asymptotic features of the time-dependent behavior. In the relaxation method of Kosloff and Tal-Ezer [7] a propagation process in imaginary time induces a relaxation of the eigenstates components with rates proportional to the eigenvalues. In the spectral method the time-energy Fourier transform of the wave function time correlation function gives an amplitude with an absolute square which peaks at the eigenvalues. [8]. The filter-diagonalization method developed in ref [9] improves the spectral method by reducing the required propagation time, permitting the resolution of small energy differences.

Direct methods involve matrix iterative procedure. The pro-

cedure consists of improving a coarse-grained solution by repeated applications of the Hamiltonian matrix or a similar matrix. The Davidson method [10] and the Lanczos algorithm [11], which fall into this category of direct methods, have been extensively used in molecular physics. These iterative methods suffer from convergence difficulties in the dense spectrum interior region where the gaps between adjacent levels may be very small. To overcome these difficulties, spectral filter methods, in which the iteration is driven by a function f(H) in place of H have been proposed within the framework of the Lanczos method. The Green function [12] exponential filter [13] and Gaussian derivative filter [14] have been successively used as filter f(H) by constructing parametrized polynominal representation of the filter functions. The Green function (E − H)^{−1} which shifts the eigenvalues of H near E to the edges of the spectrum seems to be particularly efficient. By drastically reducing the number of iteration steps in the lanczos procedure it simultaneously reduces the memory storage requirements and facilitates eigenvector calculations [15]

TheBloch wave operator theory [16] is a third direct method

which has been applied to the study of photodissociation processes using iterative treatments [17]. An advantage of the wave operator formulation is that it proposes a self-consistent strategy to reduce the dimension of the dynamical vector space by selecting the relevant active space for a given dynamical process [18] and by describing the corresponding target space. Unfortunately the various iterative solutions which have been proposed [19, 20] suffer from severe convergence difficulties. The aim of the present study is to show that significant improvements can be obtained by simply introducing into the iterative process an initial wave operator obtained by propagating the perturbed quasi-adiabatic eigenvector within the framework of time-dependent wave operator theory. We thus propose a two step procedure; the first step, related to an indirect method involves a time propagation, while the second step uses the results of the first as the starting point in a direct iterative method. Section II analyzes in detail the merits of the time-dependent wave operator formalism in performing quasiadiabatic propagations and derives some finite order differencing schemes. The performance of the method is analyzed in section III through the solution of a model system of two coupled oscillators. Other topics for investigation, such as the use of partial adiabatic representations and the introduction of spectral filters into the wave operator iterations are finally presented with the conclusions in section IV.

2 THE QUASIADIABATIC EVOLUTION OF THE WAVE OPERATOR

The procedure used to derive an eigenvector, or a group of eigenvectors within the framework of the Bloch theory has been presented in detail in several papers [21, 17]. The unperturbed vectors which are selected span an active space S_{o} with the projection operator P_{o} and simultaneously define a complementary space S_{o}^{+} with projection operator Q_{o}. The wave operator Ω = P_{o} + X, where X ≡ Q_{o}ΩP_{o}, called the reduced wave operator, is related to a non-unitary and non-singular transformation (T = 1 + X) which possesses a trivial inverse (T^{−1} = 1 − X) and which renders zero the couplings which connect S_{o} to S_{o}^{+}.

Q_{o}(1 − X)H(1 + X)P_{o} = 0  (1)

The iterative integration of the non-linear equation (1) leads to the solution of the eigenproblem. The diagonalization of H_{eff} = P_{o}HΩ inside the space S_{o} leads to the group of exact eigenvalues relating to this space.

(λ^{o})^{−1}(P_{o}HΩ)λ^{o} = P_{o}EP_{o}  (2)

The wave operator applied to the projected eigenvectors transforms them back into the exact eigenvectors

| λ_{i}⟩ = Ω | λ_{i}^{o}⟩  (3)

The active space S_{o} in which the operator Ω is constructed is itself defined by a wave operator procedure. This procedure has been explained in detail previously [17] and has been tested successfully in the study of the overtone relaxation of benzene [18] and of the CD_{3}H molecule [22]. The active space which is thus constructed and which collects together the states which are strongly coupled to the initial state undergoes a specific nonperturbative treatment: the effective Hamiltonian which corresponds to this space is stored and diagonalized at the end of the procedure.

The wave operator’s action (eq.1) is thus devoted to handling

the less strong couplings between the two spaces S_{o} and S_{0}^{+}. Nevertheless, the iterative solutions which have been proposed to solve eq.(1), such as the Recursive Distorted Wave Approximation (RDWA) [19] and the Single Cycle Method (SCM) [20] , fail when the initial state is in a dense part of the spectrum. This defect is easily understable, since these iterations make use of the ratios V_{αi}^{(n)}/(E_{i}^{(n)} −E_{α}^{(n)}) where i refers to states of the active space, α to states of the complementary space S_{o}^{+} and (n) to the order of iteration. In the dense part of the spectrum, accidental degeneracies occur during the iteration whatever choice is made for the boundaries between the two spaces. Nevertheless one can overcome these convergence difficulties by noting that the iterative series R.D.W.A. or S.C.M. procedures accept as initial input any arbitrary trial value X^{n=0}, so that the choice of a trial operator sufficiently near the required solution can make troublesome divergences disappear.

2.1 The wave operator as the adiabatic limit of the timedependent one

The time-dependent wave operator gives an efficient way to derive a coarse grained trial operator X^{n=0}. The concept of this operator was introduced to provide a simplified description scheme for the ralaxation which issues from space S_{o} and is driven by a time-dependent Hamiltonian H(t) [23]. Using notation U for the time evolution operator, we have

U(t, t_{o}; H)P_{o} = (P_{o} + X(t, t_{o}))U(t, t_{o}; H_{eff})  (4)

where

H_{eff}(t) = P_{o}H(t)(P_{o} + X(t, t_{o}))  (5)

In this description the evolution inside the space S_{o} is produced by an effective Hamiltonian (eq.5) and the transition from S_{o} to the complementary space is produced by the reduced wave operator X. The associated equation of motion for this reduced wave operator is :

ih̄\frac{∂X(t, t_{o})}{∂t} = Q_{o}(1 − X(t, t_{o}))H(t)(1 + X(t, t_{o}))P_{o}  (6)

Let use assume that the operator H is generated from the unperturbed operator H_{o} by adiabatically switching on the perturbations between t = −∞ and t = 0; the time derivative on the left hand side of eq.(6) then tends to zero and the propagation equation reduces to:

Q_{o}(1 − X(0, −∞))H(1 + X(0, −∞))P_{o} = 0  (7)

which is identical to the Bloch equation (1). The conventional Bloch operator thus appears as the adiabatic limit of the time-dependent wave operator. This offers the possibility of deriving an optimal trial operator X^{n=0} by etablishing the couplings almost adiabatically over a finite, but sufficiently large switching time T. The price to be paid is that of having to propagate with the non-linear equation (6) but fortunately this equation possesses specific features, related to its capacity to propagate quasi-adiabatic solutions, which compensate for this difficulty.

The advantages of the proposed approach are evident when one

considers the purely adiabatic evolution of a wave function issuing at t = −∞ from an eigenstate | i⟩ of the unperturbed Hamiltonian, i.e. H(t = −∞) = H_{o} and H(t = 0) = H. In the adiabatic limit the wavefunction is a product of two terms

Ψ(t) =| λ_{i}(t)⟩ exp(∫_{−∞}^{t} −iE_{i}(t^{,})/h̄dt^{,})  (8)

where | λ_{i}(t)⟩ is an instantaneous eigenvector, the components of which are slowly varying functions of time if the initial state energy E_{i}(−∞) is well separated from the other unperturbed energies, so that we can set

\frac{∂}{∂t} | λ_{i}(t)⟩ ∼ 0  (9)

The second factor in eqn (8) is a purely phase term involving the integral of the instantaneous eigenvalues.

The time behaviours of the conventional Schrödinger equation

and of the wave operator evolution equation (6) are consequently clearly different. The Schrödinger equation directly involves a finite time derivative.

ih̄\frac{∂Ψ(t)}{∂t} ≃ E_{i}(t)Ψ(t)  (10)

On the contrary the rapid phase term disappears in the expression for the wave operator constructed on the one-dimensional space S_{o} spanned by the vector | i⟩. One then has the simple result

Ω(t) = Ψ(t)/⟨i | Ψ(t)⟩
= | λ_{i}(t)⟩/⟨i | λ_{i}(t)⟩  (11)

so that:

\frac{∂Ω(t)}{∂t} ≃ 0  (12)

The smallest of the time derivative makes the time-integration easier. This feature persists when degenerate active spaces are used. These active spaces are formed in a way which supresses the largest coupling between the initial state and the surrounding bath. When for example the initial state | i⟩ is in a dense interior part of the spectrum, a convenient choice is to collect into the active space those states which are in near resonance with it. As previously noted, this artifact does not make all the divergences of the iterative process disappear in the stationary treatment, because other near resonances can arise at the boundary between S_{o} and S_{o}^{+}. This effect is less drastic in the quasi adiabatic propagation process and generally induces only some relatively small deviations from the purely stationary case. One can understand this by considering an iterative solution of eq. (6) [17]

X_{αj}^{N} = exp[\frac{i}{h̄} ∫_{−∞}^{t}(E_{j} − E_{α})dt^{,}]

× \frac{1}{ih̄} ∫_{−∞}^{t}( ∑_{n=1}^{N−1} H_{αj}^{n}(t_{1}) exp[\frac{i}{h̄} ∫_{−∞}^{t_{1}}(E_{α} − E_{j})dt_{2}])  (13)

with

H^{n} = (1 − X^{n})H(1 + X^{n}) − ih̄\frac{∂}{∂t}X^{n}  (14)

By introducing in equation (13) the approximate expression

∫_{−∞}^{t} H_{αj}(t_{1}) exp[\frac{i}{h̄} ∫_{−∞}^{t_{1}}(E_{α} − E_{j})dt_{2}]

≃ exp[i(\frac{1}{h̄} ∫_{−∞}^{t}(E_{α} − E_{j})dt_{1} + π/2)] × \frac{h̄H_{αj}(t)}{E_{α} − E_{j}}  (15)

which is valid when the time function H_{αj}(t) is effectively constant on the time scale (E_{α} − E_{j})/h̄, one obtains in the limit N → ∞ the stationary Bloch operator, with the components :

X_{αj} = ∑_{n} \frac{H_{αj}^{n}}{E_{j} − E_{α}}  (16)

in place of the time-dependent ones of (eq.13)

Thus the quasiadiabatic propagation follows the eigenvectors

continuously as long as the equation (15) is satisfied. This equation becomes incorrect when the two energies E_{α} and E_{j} are in near resonance at the border between the two spaces. However, this defect does not produce a divergence of the process such as appears in the stationary case (eq. 16) but rather induces a mixing of the instantaneous eigenvectors | λ_{α}⟩ and | λ_{j}⟩. This mixing has a small effect on a calculated eigenvector | λ_{i}⟩ situated in the middle of the active space, provided that the spectral width of this space is large enough.

2.2 The integration scheme for the time-dependent wave operator

The basic equation (6) is a non-linear equation. This means that efficient propagation algorithms such as the Chebyshev scheme [24] or the Short Iterative Lanczos propagation [25], which are effective for the linear Schrödinger equation, cannot be used. Nevertheless one can show that the use of a finite-order differencing scheme leads to expressions which, although more complicated that the equivalent expressions for the Schrödinger equation, do not require more C.P. U. integration time and can require even less if one takes into account the relative convergence rates. The second order and 4^{th} order differencing schemes are derived from the expansions

X(t + dt) = X(t − dt) + 2dtX′(t) + 0(dt^{3})
X(t + dt) = X(t − dt) + 2dtX′(t) + \frac{(dt)^{3}}{3} X′′′(t) + 0(dt^{5})  (17)

By using eq.(6) and neglecting the time derivatives of the H matrix, one obtains, after some algebraic manipulations and the renewed introduction of the wave operator Ω(t) = P_{o} + X(t), the result

X′(t) = (1 − Ω(t))C^{(1)}(t)
X′′′(t) = (1 − Ω(t)[6(C^{(1)}(t))^{3} − 3C^{(1)}(t)C^{(2)}(t) −
− 3C^{(2)}(t)C^{(1)}(t) + C^{(3)}(t)]  (18)

with

C^{(n)}(t) = H^{n}(t)Ω(t)  (19)

These propagation scheme are unitary, conserve norm and energy and, by preserving time reversal symmetry, are conditionally stable. In eqts (18) and (19), H is a full N ×N matrix Ω a rectangular N ×N_{o} matrix, where N_{o} is the dimension of the active space; the latter matrix reduces to a vector for nondegenerate actives spaces. As N_{o} is very small compared to N, the majority of the C. P. U. time used in the propagation scheme (17-19) is devoted to the construction of the rectangular matrices C^{(n)} n = 1, 2, 3. When only one vector is propagated, i.e. N_{o} = 1, this time requirement is identical to that of the finite differencing scheme for the propagation of the Schrödinger equation.

3 COMPUTATIONAL RESULTS

3.1 A two coupled oscillator model

In order to test our quasidiabatic method a simple model Hamiltonian was developed, describing a system of two harmonic oscillators in near resonance, with linear and quadratic couplings. The Hamiltonian in creation and annihilation operator form is :

H = (a_{1}^{+}a_{1} + \frac{1}{2})hν_{1} + (a_{2}^{+}a_{2} + \frac{1}{2})hν_{2}
+ λ[(a_{1}^{+} + a_{1}) + (a_{2}^{+} + a_{2}) + (a_{1}^{+} + a_{1})(a_{2}^{+} + a_{2})]  (20)

with hν_{1} = 1., hν_{2} = 1.05 and 0 ≤ λ ≤ 0.6 The unperturbed basis set is constituted of successive quasi-degenerate spaces associated with a given total number of quanta n_{1} + n_{2}. The 20 first spaces, incorporating 210 states, have been taken to constitute the representation space for the present study. In this representation we have tested our double step procedure: a quasi-adiabatic propagation leads to a trial wave operator and this zeroth-order wave operator is subsequently introduced in the iterative S. C. M. procedure to reproduce the exact wave operator. The exact eigenvectors are finally derived by using equation (3) after diagonalization of the effective Hamiltonian.

3.2 Analysis of the results

The lower part of the spectrum corresponding to the 36 first eigenvalues is presented in fig.1 as e function of the coupling parameter λ. This figure reveals that the range λ ≥ 0.3 corresponds to a strong coupling regime with many crossings and avoided crossings and thus provides a severe test for our algorithm. In the range 0. ≤ λ ≤ 0.4 the quasi-adiabatic propagation was performed by representing λ as a linearly increasing function of time on the interval [t = 0, t = N/ν_{1}] and by simply selecting as active space the quasi-degenerate space with a fixed total number of quanta equal to that for the initial state. The full spectrum corresponding to the 4 first quasi-degenerate spaces was reproduced successfully. For example, fig.2 presents the eigenvalues of the subspace number 4, with 3 quanta, versus the coupling parameter. This figure reveals one of the advantages of our procedure. In the case of parametrized Hamiltonians, it is possible to follow the eigenvectors continuously as the parameters vary and, for each new value, to establish a one to one correspondance between the zeroth-order and perturbed eigenvectors. Moreover, the construction of the full spectral map as a function of λ requires a single quasi-adiabatic propagation. Fig.3 which represents the components of the eigenvectors for λ = 0.4, exhibits a large dispersion of the results and thus confirms the strong coupling nature of the regime.

An interesting point observed in this study is the small value of

the switching time interval which is necessary to obtain a satisfactory trial wave operator. The deviations resulting from using a finite switching time interval, rather than taking the purely adiabatic case are measured by the factors.

log | ⟨λ_{j}^{exact} | λ_{i}^{app}⟩ − δ_{i,j} |  (21)

where λ^{exact} are the correct eigenvectors and λ^{app} are those obtained by the quasiadiabatic procedure. These factors are presented in figs.4 and 5 for the quasidegenerate subspace spanned by the 4 vectors with 3 quanta (3, 0; 2, 1; 1, 2; 0, 3) and for a coupling parameter λ = 0.4. The time propagation was made by using the second order differencing scheme. In the case of fig.4 the switching time interval is very small and corresponds to three classical vibrational oscillations (N = 3); it corresponds to thirty oscillations (N = 30) in the case of fig.5. As expected, one obtains a decrease of the deviation by increasing the switching time. Nevertheless the statisfactory result observed with a short switching time in fig.4 is quite remarkable; this must result from the fact that the states which are the most sensitive to the finite switching time effect do not contribute to a mixing with the initial state, since they are included in the active space. As a consequence the time propagation, which only requires 250 steps, is very fast. The use of the second order differencing scheme generates some propagation errors in this case but those remain small compared to those due to the finite switching time.

For the strongest anharmonic coupling (λ > 0.4) the S. C.

M. iteration does not converge. The analysis of fig.2 shows that our simple procedure of active space selection is not sufficient in this case, since the shifts exhibited by the levels of the quasidegenerate subspace are much greater that the distance between this subspace and the two nearest ones. In this case, convergence is obtained if adaptable and larger active spaces are seleced by the wave operator procedure.

4 SUMMARY AND DISCUSSION

The computation of small groups of interior eigenstates in a strong coupling situation is accomplished using a mixture of time-dependent and stationary wave operator approaches. The introduction of a preparatory step involving a quasiadiabatic propagation constitutes a signifiant improvement, since the direct application of the S.C.M. algorithm diverges under the same conditions. In the case of qantum relaxation processes issuing from an initial discrete state, this procedure represents an efficient tool to select the active space containing vectors which contribute significantly to the relaxation and also to calculate, in a self consistent manner, the corresponding target space. The ability of our procedure to continuously follow the vectors of this target space makes it suitable for the variational study of systems with highly parametrized Hamiltonian matrices.

The adaptability of the double-step procedure renders it ca-

pable of combination with several other methods which share the same goal. Two such methods are currently under investigation in our laboratory and have given encouraging preliminary results

The first method involves the enhancement of the convergence

by converting to a partial adiabatic representation. By using the wave-operator selection procedure it is possible to define, along with the active space a larger intermediate space which includes the latter and which collects together the states strongly coupled to the initial state. The adiabatic representation is then obtained by diagonalizing within this intermediate space. The subsequent quasiadiabatic propagation is performed by establishing adiabatically the couplings in the outer-space, the exact H-matrix projected in the intermediate space being present at the start of the propagation

A second promising method consists of using filters f(H) in

place of H in the quasi-adiabatic propagation and in the iterative treatment of the Bloch equation. The use of a filter such as the Green function (H − E)^{−1} is quite consistent with our solution procedure, since both the propagation and the iteration are based on a standard operation, namely the formation of the product of the H matrix with a vector.

Figure captions Fig1: Lower part of the spectrum of the two coupled harmonic oscillators versus the perturbation parameter λ (see. eq. 20). Fig2: The 4 eigenvalues of the quasidegenerate subspace with 3 quanta, versus the perturbation parameter λ. The eigenvalues are calculated by the double step procedure: quasiadiabatic propagation of the time-dependent wave operator and iterative S.C.M. treatment of the stationary Bloch equation. Fig3: The components ⟨j | λ_{i}⟩ of the four eigenvectors (i = 1 to 4)(∘, +, ×, ∗) of the quasidegenerate subspace with 3 quanta. The components obtained for λ = 0.4 are shown versus the index j of the harmonic basis vectors sorted in the order of increasing energy values. Fig4: The nonadiabatic factor (eq. 21) for the four states (i = 1 to 4)(∘, +, ×, ∗) of the quasidegenerate subspace with 3 quanta. The factors are shown versus the index j of the harmonic basis vectors sorted in the order of increasing energy values. The switching time of the perturbation is equal to 3 vibrational periods. Fig5: The same as fig.4 but with a switching time equal to 30 vibrational periods.

References

[1] Junker, B.R. 1982. Advan. At. Mol. Phys. 18 : 207

[2] Reinhardt, W.P. 1982. Ann. Rev Phys. Chem. 33 : 223

[3] Kukulin, V. I. Krasnopol’sky, V. M., Horacek, J. 1989. Theory of resonances. Dordrecht : KLuwer Academic Publishers

[4] Baslev, E., Combes, J. M. 1971. Commun. Math. Phys. 22 : 280; Aguilar, J., Combes, J.M. 1971. Commun. Math. Phys. 22 : 265; Simon, B. Commun. Math. Phys. 27 : 1

[5] Moiseyev, N., Hirschfelder, J. O. 1988. J. Chem. Phys. 88 : 1063

[6] Jolicard, G., Austin, E. J. 1985. Chem. Phys. Letters 21 : 103; Jolicard, G., Austin, E. J. 1986. Chem. Phys. 103 : 295

[7] Kosloff, R., Tal-Ezer, H. 1986. Chem. Phys. Letters. 127 : 223

[8] Feit, M. D., Fleck, J. A., and Steiger, A. 1982. J. Comp. Phys. 47 : 412

[9] Neuhauser, D. 1990. J. Chem. Phys. 93 : 2611 ; 1991 J. Chem. Phys. 95 : 4927

[10] Davidson, E. R. 1989. Comp. Phys. Comm. 53 : 49

[11] Lanczos, C. 1950. J. Res. Natl. Bur. Stand. 45 : 255

[12] Ericsson, T., Ruhe A. 1980. Math. Comp. 35 : 1251

[13] Webster, F., Rossky, P. J. , Friesner, R. A. 1991. Comp. Phys. Comm. 63 : 494

[14] Kono, H. 1993. Chem. Phys. Letters 214 : 137

[15] Wyatt, R. E. Phys. Rev. E. unpublished

[16] Bloch, C. 1958. Nucl. Phys. 6 : 329

[17] Jolicard, G., Périé, J. 1993. Advances in multi-photon processes and spectroscopy Vol 8, p. 1-126 Singapore : World Scientific

[18] Wyatt, R. E., Iung, C., Leforestier, C. 1992. J. Chem. Phys. 97 : 3458, ; ibid 97 : 3477

[19] Jolicard, G. 1987. Chem. Phys. 115 : 57

[20] Killingbeck, J. P., Jolicard, G. 1992. J. Phys. A Math. Gen. 25 : 6455 [21] Barthelat, J. C. Durand, Ph. 1992. In Computational Chemistry. Struc-

tures -interactions and Reactivity Part A : p.26, ed S. Fraga Amsterdam : Elsevier

[22] Iung, C., Leforestier, C. J. Chem. Phys. unpublished

[23] Jolicard, G. 1989. J. Chem. Phys. 90 : 2320 ; Jolicard, G., Killingbeck J. P. 1994 Domain-based parallelism and problem solving decomposition methods in computational science ans engineering S. I. A. M. Book, in press

[24] Tal-Ezer, H., Kosloff, R. 1984. J. Chem. Phys. 81 : 3967

[25] Park, T. J., Light, J. C. 1986. J. Chem. Phys. 85 : 5870