• Keine Ergebnisse gefunden

C omputational M ethods

4.2 Methods of nuclear quantum dynamics

expressed as sum of an effective Hamiltonian, whose eigenstates are the CASSCF wavefunctions, and a perturbation operator. The wavefunction and the electronic energy are expanded up to the first and second perturbative order, respectively. Only single and double electronic excitations from the CASSCF wavefunction are necessary. Rayleigh-Schr ¨odinger perturbation theory is used to derive the combination coefficients and the second-order corrected energy. In this work, CASPT2 calculations were used to determine the potential energy surfaces of the ground and theπσstates of pyrrole.

4.2 Methods of nuclear quantum dynamics

The quantum mechanical study of the nuclear motion requires the solution of the time-dependent Schr ¨odinger equation (2.24) for the nuclei,

i¯h|Φ,ti

∂t = H |Φ,ti , with |Φ, 0i=|Φ0i . (2.24)

As discussed in Section (3.2), the wavefunction Φ for the atom–fragment photodissociation is labelled with the total angular momentum quantum numbers (J,M)and depends, in general, on the Jacobi coordinates (R,θ,φ), the Euler angles Ω, and the fragment internal coordinatesQ. In this Section, the coordinates are generically indicated as xκ, κ=1, ...,f.

With modern computer architectures, the TDSE can be solved exactly for low-dimensional problems (up to five or six vibrational degrees of freedom), using standard polynomial methods, such as Chebyshev propagation.51In order to propagate the wavefunction for high-dimensional systems, several accurate techniques have been developed,37among which the most popular is the multi-configurational time-dependent Hartree method (MCTDH).82,83

4.2.1 Grid representation and Chebyshev Propagation

For low-dimensional systems, such as a triatomic molecule, the wavefunction can be represented on a multi-dimensional grid, typically using a discrete variable representation (DVR).84If a direct tensor product grid is used, the wavefunction has the form

Φ(x1, ...,xf) =

N1

j

1=1

...

Nf j

f=1

Cj1...jf(t)θ(j1)

1 (x1...·θ(jff)(xf), (4.1)

whereθ(jκ)

κ is the DVR function for the coordinatexκ, localized at the jκ-th point on the grid for xκ. The integersNκare the grid dimensions and the complex numbersCj1...jf(t)are time-dependent coefficients associated with the multi-dimensional grid points.

Using the representation of Equation (4.1), the TDSE can be written in the form i¯h∂CJ

∂t =

K

HJLCL, (4.2)

where J = (j1...jf) is a multi-index and HJL = Dθ(j1)

1 ·...·θ(jff)Hθl(1)

1 ·...·θ(lff)E is the matrix-representation of the Hamiltonian in the (multi-dimensional) DVR basis.

Equation (4.2) can be solved using efficient polynomial methods: With H = {HJL} and C={CJ}, Equation (4.2) is solved as

i¯h∂C

∂t =HC =⇒ C(t) =exp

i

¯ hHt

C(0). (4.3)

The propagator is represented as exp −h¯iHt

=nanPn(H), wherePn(H)is a polynomial of the Hamiltonian matrix.

In Chebyshev propagation, the Hamiltonian matrix is first rescaled to the domain[−1, 1]: H −→ Hnorm= 2H1(Emax+Emin)

EmaxEmin

, (4.4)

whereEminand Emax are estimations for the minimum and maximum eigenvalues ofH. Then, the coefficient vector is calculated as51,85

C(t) =eih(Emax+Emin)t

Niter

n

=0

(−i)n(2−δn0)Jn

EmaxEmin

2¯h t

Cn , (4.5)

where the coefficients are expressed in terms of the Bessel functions Jn(α), and the vectorsCnare generated by the recurrence relation:

C0 =C(t=0) C1 =HnormC0

Cn+1 =2HnormCnCn1 .

(4.6)

Equation (4.5) shows that the time-dependence of the wavefunction is given by the Bessel function coefficients and the spatial dependence is in theCn vectors. Autocorrelation or cross-correlation functions are calculated by storing the overlaps of theCn vectors either withC0 or with the vectors of interest (associated, for example, with fragment eigenstates); such time-independent overlaps are then weighted with the time-dependent Jn coefficients whose Fourier transform is obtained analytically.

4.2 m e t h o d s o f n u c l e a r q ua n t u m d y na m i c s 35

The computational effort of this numerically exact method is mainly given by cost of the HCn matrix-vector multiplication. Assuming, for simplicity, N1 = ... = Nf = N, and that the kinetic energy operator is given as a sum of single-coordinate terms and the potential operator is diagonal on a DVR grid, the computational cost is proportional to f Nf+1.86 This exponential scaling, common to all polynomial integrators, limits the use of the method to systems with few degrees of freedom.

4.2.2 The multi-configurational time-dependent Hartree method (MCTDH)

The MCTDH method83,82is based on the representation of the wavefunction in terms of configura-tions. A set of time-dependent single-particle functions (SPFs) is used for each degree of freedom, and the configurations are built as Hartree products of SPFs. SPFs are typically represented on a DVR grid and the exact description of Section 4.2.1 is recovered when the number of SPF approaches the size of the DVR basis.

The time-dependent wavefunction is given by the followingansatz:

Φ(x1, ...,xf,t) = coefficients. In the MCTDH wavefunction, not only the coefficients, but also the ‘basis functions’

φ(jκ)

κ depend on time.

The equations of motion are derived from the time-dependent variational principle,

The variationsδΦmust be performed with respect to the coefficients, δΦ where the single-hole functions Φe(Jκ) have been defined as combinations of Hartree products between (f−1)SPFs.

Inserting Equations (4.9) and (4.10) into Equation (4.8), the MCTDH working equations are obtained after some cumbersome rearrangements:

i¯h∂CJ

Equation (4.12) requires the evaluation of the projection operator P(κ)=

and of the product between the inverse of the density matrix ρ(jlκ) = DΦe(jκ)(lκ)E

The computational cost for the MCTDH algorithm is the sum of two parts.83One contribution is due to the action of the mean field matricesD

H(k)E

jl on the SPF functions. It grows linearly with the number of degrees of freedom f and with the average numbernof SPFs for each degree of freedom and quadratically with the average dimension N of the one-particle DVR grid. The total effort scales as∼n f N2. The second contribution is due to the cost of building the mean field matrices and grows as∼ f2nf+1. The total computational cost is therefore

computational cost≈c1f nN2+c2f2nf+1 , (4.13) with coefficients of proportionality,c1andc2. The effort of the standard method of Section 4.2.1is proportional to f Nf+1. It follows that, for large values ofnand f, the gain factor of the MCTDH

In order to remove redundancies between SPFs and coefficients, the constraints D

4.2 m e t h o d s o f n u c l e a r q ua n t u m d y na m i c s 37

scheme is proportional to f1(N/n)f+1: For a high number f of degrees of freedom and for a large mean contraction efficiency(N/n), the MCTDH method is much more efficient than the standard wave packet propagation method.

5