A Practical Diabatisation Scheme for Use with the Direct-Dynamics

Article pubs.acs.org/JPCA

A Practical Diabatisation Scheme for Use with the Direct-Dynamics Variational Multi-Configuration Gaussian Method Gareth W. Richings and Graham A. Worth* School of Chemistry, University of Birmingham, Edgbaston, Birmingham B15 2TT, U.K. ABSTRACT: A method for diabatising multiple electronic states on-the-fly within the direct dynamics variational multiconfiguration Gaussian method for calculating quantum nuclear dynamics is presented. The method is based upon the propagation of the adiabatic−diabatic transformation matrix along the paths followed by the Gaussian basis functions that constitute the nuclear wave function, by use of a well-known differential equation relating the matrix and the nonadiabatic vector coupling terms between the electronic states. The implementation of the method is described, and test calculations are presented using the ground and firstexcited states of the butatriene cation as an example, allowing comparison to the earlier regularisation diabatisation scheme as well as to full nuclear dynamics on a precomputed potential energy surface. The new scheme is termed propagation diabatisation.

■

INTRODUCTION The theoretical study of nuclear dynamics in molecules is of fundamental importance in the understanding of scattering processes, both reactive and nonreactive, spectroscopy, dissociation, etc. The feature common to any such study is a representation of the nuclear wave function moving in a potential described by the electrons. Such simulations of the dynamics of nuclei can be performed using classical mechanics, but where the motion of the wave function occurs on multiple, generally interacting, potential energy surfaces (PESs), which represent distinct electronic states, it is necessary to include the effects of quantum mechanics,1−12 either rigorously or in an ad hoc manner, to account for the possibility of transfer of parts of the wave function between the surfaces. The traditional method for performing nuclear dynamics calculations is to generate the PESs prior to performing the actual nuclear dynamics. The methods for producing such surfaces range from using model potentials, either theoretically or empirically generated, through to the use of potential energies generated at assorted molecular geometries using an appropriate electronic structure method (ab initio or semiempirical). PESs can be generated from such sets of electronic energies by fitting analytic functions of a predetermined form to the data set, as exemplified by the vibronic coupling Hamiltonian method,5,13 which uses a harmonic approximation to the potential as the starting point before adding extra terms to produce a more accurate fit, or by spline-fitting between the points. The nuclear dynamics can be very sensitive to the actual PES used, so for example, when trying to reproduce and explain an experimentally derived spectrum, the precision with which the PES is calculated can have a great bearing on the success of © XXXX American Chemical Society

such an endeavor. This means that high-level ab initio calculations, possibly involving many electronic states, may be needed. If many nuclear degrees of freedom are needed for an accurate picture, this will require many calculations. The fitting of the surfaces then needs to be very accurate, meaning that complicated functions with many parameters must be used to represent the PES.14−16 This can be a very time-consuming process, and must be completed before the dynamics calculation can proceed. To alleviate this bottleneck a lot of recent work has focused on so-called direct dynamics methods where the PESs are created on-the-fly, that is, as the nuclear dynamics proceeds (recent reviews can be found in refs 17 and 18). This has the advantage that the formation of the PES is governed by the dynamics: the motion of the nuclei determining where, in configuration space, more information about the PES is needed. Such methods include those based on the trajectory surface hopping (TSH) method,19−26 where the wave function is represented by a swarm of classical trajectories moving over the PES. These calculations can be performed on fitted surfaces, but in the on-the-fly version electronic structure calculations are performed along the trajectories.27−36 In both cases transfer between states is determined by probabilities based around the Landau−Zener theory of curve crossing, the most popular method for doing so being Tully’s fewest switches algorithm.20 Special Issue: Dynamics of Molecular Collisions XXV: Fifty Years of Chemical Reaction Dynamics Received: August 14, 2015 Revised: September 29, 2015

A

DOI: 10.1021/acs.jpca.5b07921 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

review,66 but we will give a brief outline of the parts necessary to develop what follows for completeness. The vMCG method is one of several ways of solving the time-dependent Schrödinger equation

However, to gain accurate results many trajectories may be needed, with the inherent increase in computational cost that entails, and there is no guarantee of convergence on the correct result. The ab initio multiple spawning method (AIMS) of ́ 37−40 (the direct dynamics version of the multiple Martinez spawning method41,42) improves on TSH by replacing the pointlike trajectories with a wave function constructed as a linear combination of time-dependent Gaussian functions. The Gaussians move classically, but their time-dependent coefficients are determined quantum mechanically. To effect crossing between PESs, nonadiabatic couplings between the states are calculated at the centers of the functions, and once they exceed some threshold, new Gaussians are spawned on the other surface at the same geometry. The method is fully quantum mechanical but due to the classical basis functions can suffer from the same convergence issues as TSH. To improve convergence a fully quantum method is needed, and one such method will be the focus of the work presented here. The method is the direct-dynamics variational multiconfiguration Gaussian method (DD-vMCG), which, like AIMS, uses a linear combination of Gaussian wavepackets (GWPs) to represent the wave function. The key difference from AIMS is that the motion of the GWPs is described fully quantum mechanically. Many studies of molecular systems have been performed using DD-vMCG including benzene,43,44 nitrosyl chloride,45 fulvene,46,47 cyanine,48 thymine,49 and formaldehyde.50,51 As will be more fully explained in later sections, the DDvMCG method relies on having a PES where all first and second derivatives are defined. This, however, is not the case at points in the molecular configuration space where electronic states become degenerate. In the adiabatic representation, where the states are ordered in increasing energy, the relevant PESs have an undefined first derivative at such points, and the nonadiabatic couplings become infinite. To alleviate this it is necessary to transform to the so-called diabatic representation where all PESs are adequately smooth at all points.12,52−54 Various ways of defining diabatic states have been proposed over the years,55−63 and the problem is not straightforward when fitting functions to a precalculated set of electronic structure energies. It is even less so when constructing the PESs on-the-fly. The aim of this work is to present a solution to this problem in the paradigm of DD-vMCG, applicable to any number of states and to any number of state degeneracies. In the next section we will describe the theory behind DDvMCG and the diabatisation schemes available to us. Subsequently we present the new diabatisation method for DD-vMCG along with details of its implementation in the Quantics package.64 The penultimate section presents test calculations on reduced-dimensionality models of the butatriene molecule, a system that was studied previously by one of us in relation to TSH and DD-vMCG. Finally, conclusions along with some ideas for future development will be presented.

iℏ

∂Ψ(R, t ) = Ĥ Ψ(R, t ) ∂t

(1)

for a wave function Ψ, which is a function of both nuclear coordinates, R, and time, t. The method uses a linear combination of multidimensional Gaussian functions (GWPs), {gj}, as an ansatz. We will be interested in dynamics on multiple PESs and will be using the so-called single-set formulation of vMCG, where a single set of GWPs is used to expand the wave function on each state, this being done by having a separate set of coefficients for each state. The ansatz is N

Ψ(R, t ) =

Ns

∑ gj(R, t ) ∑ A(j s)(t )|s⟩ j=1

s=1

(2)

{A(s) j }

where is the set of coefficients for the state s. Crucially, both the coefficients and the GWPs are time-dependent, meaning that the weight of each GWP can vary over time and also that the positions and momenta of the GWPs can evolve in molecular phase space. This means that the basis can move with the evolving wave function, and be kept optimally small. Insertion of this ansatz into the Dirac−Frenkel variational principle ⟨δ Ψ|H − iℏ

∂ |Ψ⟩ = 0 ∂t

(3)

yields two sets of equations-of-motion (EOMs), one for the set of coefficients and one for the parameters of the GWPs that allow them to evolve variationally (determining position, momentum, width, and phase). See ref 66 for details. Given an appropriate representation of the PESs, solution of the EOMs using appropriate initial conditions and standard numerical integrators then allows the nuclear dynamics on those surfaces to be followed. The vMCG method is of general utility for nuclear dynamics and can be used on standard fitted PESs, but it is of most use when combined with the concept of direct dynamics, to give the DD-vMCG method. In this case the PESs, instead of being determined prior to performing the dynamics, are created onthe-fly. To do this, an electronic structure program is used to calculate energies of the system in question at molecular geometries determined by the positions of the GWPs as they move in configuration space. As the GWPs move, energies and the associated gradients and Hessians are calculated at the centers of each of them, and this information is stored in a database. The gradients and Hessians of the electronic energies are needed for two reasons: first to allow the evaluation of potential energy matrix elements in the vMCG EOMs by use of a second-order expansion around the GWP center (local harmonic approximation). Second, the expansion of the energy around the geometries in the database means that it is not necessary to recalculate energies, gradients, and Hessians at every time step during the propagation of the nuclear wavepacket. New database points are only calculated once the molecular geometry, represented by the GWP center, differs by some predetermined amount from all database geometries. While the GWP is close enough to the database points, energies, gradients, and Hessians are interpolated from the

■

COMPUTATIONAL DETAILS Direct-Dynamics Variational Multi-Configuration Gaussian. As explained in the Introduction, we aim to present a diabatisation scheme for use with the DD-vMCG method65 within the Quantics package. Full details of this method and the underlying vMCG dynamics approach can be found in a recent B


Article


positions of the centers of 25 GWPs for a two-dimensional model of butatriene (which will be more fully described later). The coordinates are normal modes that belong to two different irreducible representations of the D2h point group. Because of the properties of the Au representation, the PES is symmetric about the x-axis of this plot and so should the wavepacket remain if it begins centered on that axis. In this example the GWPs labeled 2 and 4 are symmetry-equivalent, and as a result both have an equivalent overlap with, for example, GWP 3. This equivalence should follow through to the inverse of the overlap and C-matrices, but it does not quite do so. The result of this is an accumulation of errors that causes the wavepacket to move to one side or other of the axis. Various inversion methods and regularisation schemes were used to try and make the EOMs behave, to no avail. It was thus decided to enforce symmetry on the wavepacket after each step in the propagation. This was done by noting symmetryequivalent GWPs at the start of the propagation. At the end of each step, the positions and momenta of the equivalent GWPs were checked and then equalized to the mean of the absolute values (the relative signs then being maintained). Also noted were those GWPs that began on the symmetry plane, for example, GWPs 1, 3, 5, 11, and 13 in Figure 1, and which must stay on the plane during the propagation. Their coordinates were constrained to be zero in the relevant direction. The breaking of symmetry is also reflected in the A-coefficients (eq 2) once the GWPs move away, so these also need to be equalized between equivalent GWPs, and it is also necessary to ensure that the GWPs constrained to the plane of symmetry have zero coefficients if the wavepacket has a node there (determined by checking the phase change across the plane). Adiabatic and Diabatic Representations of Potential Energy Surfaces. The PESs generated by performing electronic structure calculations at various molecular geometries are so-called adiabatic states, {|ψi⟩}. These surfaces connect points with the energies in ascending order, for example, the third adiabatic state is the surface connecting the third lowest energy at all geometries. If we use these states as a basis for expanding the full wave function, with coefficients provided by the nuclear wave functions {χi}, then

values present in the database. The method by which this is done is based on the GROW philosophy of Collins and coworkers67,68 as adapted by Frankcombe et al. for vMCG functions.69 In effect the interpolated values are generated by a Shepard interpolation as V (R) =

∑i ∈ database Ti |R − R i|−2p ∑j ∈ database |R − R j|−2p

(4)

where Ti is the Taylor series expansion of the electronic energy, gradient, or Hessian (second-, first-, and zeroth-order, respectively) centered at the database point Ri, and p is an integer, for which we choose p = 2. As we are interested in dynamics on more than one PES we will need to calculate excited-state energies, etc. In our example calculations, presented later, this will be done using the complete active-space self-consistent field method (CASSCF). As is well-known, the choice of active space is vital for the success of such calculations. The active space is chosen at the molecular geometry where the dynamics begins, and when the electronic structure calculation is performed, the resultant molecular orbital (MO) coefficients are stored in the database. When a new calculation must be performed, the coefficients from the database point that has the geometry most similar to the new point are read and used to give the initial guess wave function. The resulting MO coefficients are then stored in the database for later use. It was noted during the initial calculations on butatriene for this work that the nuclear wavepacket invariably broke the point group symmetry of the molecule. Where the wavepacket should have had a plane of symmetry, this was lost very early on in the propagation (typically within the first femtosecond), with the deviation from the plane being of the order of 0.7 (in unitless, mass-frequency scaled coordinates) after 5 fs. The reason for this was traced to numerical instabilities in the inversion of the matrix of GWP overlaps (necessary due to the nonorthogonality of the GWP basis) and of the so-called Cmatrix in the EOMs of the GWP parameters (see ref 66). To symmetrize the wavepacket, the following procedure was implemented. Consider Figure 1 where we have the initial

|Ψ(R, r)⟩ =

∑ |χi (R)⟩|ψi(r; R)⟩

(5)

i

where R is the vector of nuclear coordinates, and r is the vector of electronic coordinates. The electronic basis is parametrically dependent on the nuclear coordinates. By insertion of this wave function into the time-dependent Schrodinger equation and projection from the left by a member of the set of electronic basis functions, ψj, we get the Schrödinger equation for the nuclear wave functions ̂ + Vjĵ )|χ ⟩ − (Tnuc j

∑ Λ̂ji|χi ⟩ = iℏ i

∂|χj ⟩ ∂t

(6)

where T̂ nuc is the nuclear kinetic energy operator, V̂ jj is the potential due to the electronic state j, and {Λ̂ji} are the nonadiabatic coupling operators between states i and j. The latter are due to the operation of the nuclear kinetic energy operator on the electronic basis functions, which is dependent on the nuclear coordinates and has the form

Figure 1. Initial positions of the centers of 25 GWPs, as determined by the algorithm implemented in the Quantics package, constituting a wavepacket centered at the origin, along the unitless mass-frequency scaled normal modes 5Au (torsion) and 14Ag (central C−C stretch) of butatriene. Labels of the positions reflect the order in which additional GWPs are added when the total number is increased.

Λ̂ji ∝ (2Fji ·∇ + Gji) C

(7) DOI: 10.1021/acs.jpca.5b07921 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article


system of two states we can define the following quantities in terms of the adiabatic energies

where the constant of proportionality is related to the inverse of the nuclear masses. Gji is the so-called scalar coupling term, and Fji is the vector coupling term, which has the form Fji =

⟨ψj|∇Ĥ |ψi ⟩ Vii − Vjj

=

ψj ∇ψi

1 A (V22 − V11A) 2

(12b)

⎛−Δ V D ⎞ D 12 ⎟ V D = ΣI + Δ·⎜⎜ ⎟/Δ D ⎝ V12 ΔD ⎠

(13)

(14)

where ΔD is the diabatic analogue of Δ. Clearly this equation becomes singular at a point of degeneracy of the adiabatic states, so the assumption is made that the adiabatic and diabatic states are equal at such a point (a valid condition as will be more fully described later). By doing so, the diabatic matrix can be expanded to first-order in terms of the gradient difference vector of the adiabatic states, ∂ = ∇(VA22−VA11), where the derivatives, here and in all that follows, are with respect to the nuclear coordinates, and the derivative coupling vector between them, F12.

(9)

⎛ ∂ F12 ⎞ ⎜⎜ ⎟⎟ ·dR V (1) = D ⎝ F12 ∂ ⎠

(15)

where dR is the nuclear displacement from the point of degeneracy. The original method, as applied to predetermined sets of adiabatic energies, uses molecular symmetry to uniquely define the elements of the matrix in eq 14,71 the diagonal elements only containing components from totally symmetric molecular modes, and the off-diagonal elements consisting of elements that transform according to the direct product of the irreducible representations of the interacting states. The method used in DD-vMCG takes the orthogonal gradient difference and derivative coupling vectors, calculated at a point of degeneracy to provide this uniqueness. The adiabatic energy gap can be approximated by Δ(1) = |∂|2 + |F12|2 , which is a positive number, and thus the diabatic states can be represented

∑ |ψj(R)⟩cji(R) (10)

V D = ΣI + Δ·

Note that all components of this relation are dependent on the molecular configuration space, that is, the states and the transformation between them are functions of molecular geometry. For clarity, we will not write this explicitly in the equations that follow. At each point, the adiabatic and diabatic energy matrices are related by a similarity transform using the matrix constituting the coefficients in eq 10, C. V D = C†V A C

Δ=

where I is the unit vector. As the transformation is unitary, the trace of the transformed matrices is invariant, and hence 1 D D Σ = 2 (V11 + V22 ). This leads to

N j=1

(12a)

⎛− 1 0 ⎞ ⎟C V D = ΣI + Δ·C†⎜ ⎝ 0 1⎠

in general. As such we get smooth surfaces through all points of degeneracy. It can be shown that the vector coupling can only be made rigorously zero for diatomics52,54 or if an infinite set of adiabatic states is available.12 As a result we must settle for quasi-diabatic states {|φi⟩}, where the couplings are minimized enough to be ignored and to still give smooth PESs. It is with these states that we will be concerned in what follows, where we will omit the prefix, quasi, for brevity. Regularisation Diabatisation. The diabatisation scheme used in DD-vMCG to date is based on Köppel’s regularisation diabatisation method.70−72 We will briefly summarize the main points, following refs 65 and 66, as a way of comparison with the new scheme presented in the next section. The set of N quasi-diabatic states {|φi⟩} is related to the set of adiabatic states {|ψi⟩} by the unitary transformation defined as |φi(R)⟩ =

1 A A (V11 + V22 ) 2

and hence

(8)

The second equality comes from consideration of the Hellman−Feynman theorem, which holds for the adiabatic basis, and from the fact that ∇⟨ψj|Ĥ |ψi⟩ = 0 everywhere in nuclear configuration space. The vector coupling term, being inversely proportional to the energy gap between the ith and jth adiabatic states, is safely ignored when that gap is large, but becomes significant when the states become close in energy. Worse though is that, when the states become degenerate, the coupling vector becomes singular. In effect, at such a point the gradient of the adiabatic states becomes undefined, which is a problem if we are trying to perform nuclear dynamics on such a surface. To alleviate these issues, it is possible to transform the adiabatic states to the so-called diabatic representation. Pure diabatic states {|ϕi⟩} are (geometry-dependent) linear combinations of the adiabatic states constructed so that the singular coupling terms vanish,12,52 the coupling being transferred to well-behaved potential like terms, that is:

⟨ϕi|Ĥ |ϕj⟩ ≠ 0

Σ=

V (1) D Δ(1)

(16)

This is the equation that gives the diabatic states in the regularisation diabatisation method as implemented here. By construction, the adiabatic energies can be retrieved by diagonalization of the diabatic energy matrix, but the linear approximation introduces errors into the nonadiabatic couplings that are only exact at the point of degeneracy. The method suffers from three major drawbacks: the first is that it is only applicable to system of two states; in general, we are interested in the dynamics of systems involving more than two. The second problem in the context of direct dynamics is that the location of the degeneracy must be found before the

(11)

where VA is the diagonal matrix of adiabatic energies, and VD is the, in general, nondiagonal matrix of diabatic energies. For a D


Article

The Journal of Physical Chemistry A dynamics can be performed. This rather goes against the spirit of a direct method, and may also be a nontrivial task. The third issue is that, in more than two dimensions, there are an infinite number of points of degeneracy. How then should the reference point be chosen? To date, the minimum energy point on a conical intersection seam has been taken as representative of the seam, a good approximation if this is the crossing point of the wavepacket. In the next section we present a diabatisation scheme based on the propagation of the adiabatic to diabatic transformation matrix,53 which does not suffer from these drawbacks: it can in principle deal with an arbitrary number of states, needs no prior information about the PESs, and can deal with any number of crossings in an even-handed way. Propagation Diabatisation. Considering eq 10 we can use the fact that the adiabatic states form an orthonormal set, to show that the elements of the transformation matrix, C, are simply overlaps of the relevant adiabatic and diabatic states cji = ⟨ψj|φi⟩

Fij =

∇cji =

k=1

k=1

⎛1 exp⎜ ⎝2

(19)

N

∑ ⟨∇ψj|ψk⟩cki k=1

(20)

Defining the nonadiabatic coupling vector as in eq 8, it is easy to see, by remembering the orthonormality of the adiabatic states, that N

∇cji ≈ − ∑ Fjk cki k=1

(21)

and in matrix form ̲ ∇C ≈ −FC

∫R

R + ΔR

F̲ ·dR)C(R)

(24)

⎞

⎛

⎞

∫ F̲ ·dR⎠C(R + ΔR) = exp⎝− 12 ∫ F̲ ·dR⎠C(R) ⎟

⎜

⎟

The exponentials are then expanded in Taylor series, and by inversion of the resultant matrix on the left-hand side, we get the final transformation matrix, C(R + ΔR). This then is the essence of the propagation diabatisation method that we propose for use with DD-vMCG. Our method follows a similar philosophy to the local diabatisation scheme of Granucci et al.75 in that the form of the diabatic states is propagated along with the nuclear wavepacket basis functions (in their case, the classical trajectories of TSH). In that method the diabatic states are defined to be those where the dynamic ∂ couplings ⟨φi| ∂t φj⟩ (which are related by the chain rule to the vector couplings) between two diabatic states vanish. This only occurs along the velocity vector of the trajectory in question, hence the local nature of the method. Propagation of the coefficients of the diabatic states in which the electronic wave function is expanded allows the determination of the adiabatic− diabatic transformation matrix. This method has been modified and applied to the SHARC variant of TSH,26 allowing propagation of the so-called diagonal representation of the PESs (where the potential matrix including arbitrary couplings such as spin−orbit couplings is diagonalized) as well as of the transformation to a diabatic representation.76,77 The localized diabatisation method thus uses propagation of the adiabatic− diabatic transformation matrix through time, whereas the propagation diabatisation method integrates through molecular configuration space. However, both methods assume that couplings to states outside the manifold under consideration can be ignored and are both dependent on the path traveled by the basis functions, features that will be more thoroughly discussed in the next subsection. Our implementation of the propagation diabatisation method will be described in the following sections. Removable and Non-Removable Couplings. As noted in the section Adiabatic and Diabatic Representations of Potential Energy Surfaces the vector coupling terms can only be

This equation holds trivially for the so-called crude adiabatic basis, where the electronic basis is independent of nuclear geometry, but this form also holds for more general bases. As we are only dealing with a quasi-diabatic representation, this will not strictly hold, but if we suppose the couplings are sufficiently small, in some sense, to be ignored then12,53 ∇cji ≈

(23)

(25)

(18)

In a strictly diabatic representation, the final term1,52 ⟨φk |∇φi⟩ = 0, ∀ i , k

Vjj − Vii

To solve this equation, we are faced with two main problems, namely, the integral of a function with unknown analytic form and the exponential of a matrix expression. The first problem is simply overcome in this proof-of-concept work by a numerical integration using the trapezium rule along a straight line between R and R + ΔR. Doing so we find a matrix of scalars, of which we must take the exponential. The second problem requires us to consider the subtlety that the propagation, as expressed in eq 24, does not guarantee that a unitary matrix at the start will yield a unitary matrix at the end.74 To maintain the unitarity of the transformation matrix we follow the method of Esry and Sadeghpour74 where they used the Cayley−Hamilton form of the propagator. In essence we rearrange eq 24 as

N

∑ ⟨ψj|φk⟩⟨φk|∇φi⟩

Dij

C(R + ΔR) = exp( −

(17)

∑ ⟨∇ψj|ψk⟩⟨ψk|φi⟩ +

Vjj − Vii

:=

The solution of this differential equation (22) defines our diabatisation scheme. The equation has the formal solution of a propagation over some short step,53 ΔR, starting at position R

As this is a scalar function of position, we can take the gradient of the matrix element, and by appreciating the fact that both the sets of adiabatic and diabatic states are individually complete, we see that N

⟨ψi|∇Ĥ |ψj⟩

(22)

where the underlining reminds us that F is a matrix of vectors. This equation was first presented by Baer,53 who investigated the necessary and sufficient conditions for the equality to have a solution, along with its transformation to an integral equation. It is this integral equation upon which we will be basing our proposed scheme. For future use, and recalling eq 8, we note that the nonadiabatic coupling vector is the ratio of the derivative coupling vector (the quantity generated by the electronic structure program, Gaussian 03,73 used here) between the states in question (i and j) and their adiabatic energy gap, that is: E


Article


removal of the infinite couplings occurs independent of which states are involved within the space being described (two or more). Another point to note regarding the removability of couplings is that, because eq 22 is not a strict equality, the integrals in eq 25 are not path-independent.54 In the next section, we will indicate the choice of path that we make in the practical implementation of the propagation diabatisation method. The path dependence will be weakest in the vicinity of any intersections between states included in the calculations because the removable couplings dominate, however many states are included in the manifold under consideration, and eq 22 will hence be an approximate equality in these regions. As such, the description of the PES around the conical intersection, the region of most importance to treat accurately, will be least affected by this potential source of inaccuracy. Away from any intersections we assume that all nonadiabatic couplings are negligible, that is, we use the Born−Oppenheimer approximation and do not propagate the transformation matrix. Implementation. To test the propagation diabatisation scheme, it was coded into a development version of the Quantics package.64 This program, which grew out of the wellknown Heidelberg MCTDH Package,78 has the capacity to run DD-vMCG calculations on multiple PESs. The Quantics package is interfaced with the Gaussian 0373 electronic structure program to provide the necessary energies, gradients, Hessians, and derivative couplings. The implementation is as follows: Initial Conditions. An initial nuclear wavepacket is defined, centered at a particular molecular geometry. This wavepacket is constructed as a linear combination of multidimensional GWPs as described in the Computational Details Section “DDvMCG”. The location of the basis functions is determined automatically, but in general it is done by building concentric shells of GWPs around the center of the wavepacket. In placing the GWPs, a balance must be struck between an accurate representation of the wavepacket and the linear dependence of the set of GWPs (which are not orthogonal to one another). Considering eq 22, we see that, by multiplying C by any unitary matrix, which is constant over the configuration space, the form and hence solution of the equation is unchanged. As a result, we are free to choose the value of C at the center of the wavepacket by a judicious choice of this constant matrix. We therefore choose Cref. = I, the unit matrix. This means that at the center of the initial wavepacket, the adiabatic and diabatic representations of the electronic states are assumed to be equal. As noted in Section DD-vMCG, electronic structure calculations are performed at the center of the GWPs. Usually, at t = 0, there will be a GWP at the center of the wavepacket, which can be used as the reference point. However, if the widths of the GWPs differ from that of the wavepacket, there will not necessarily be a GWP at the center of the wavepacket (so as to most accurately represent the chosen form of the wavepacket). In this case an additional, reference database point is created at the center of the wavepacket, before subsequent database points are generated at the centers of the GWPs. Transformation of Energies and Gradients. Having calculated the adiabatic energies, gradients, Hessians, and derivative couplings at the central point (for each adiabatic state), we can transform them to the diabatic representation that is used for the nuclear dynamics. Obviously this is trivial for the case where C = I, but the principle holds for all C. The

completely removed (for a nonlinear molecule) if we have a complete set of adiabatic states. As this is obviously an impossible scenario, it is worth considering the situations in which eq 19 holds, at least to a good degree of approximation, meaning that eq 22 can be solved numerically and give reasonable results.54 The conditions under which eq 22 has an exact solution have been thoroughly studied previously.12,53 Earlier work54 has divided the vector coupling terms in a finite basis of states into removable and nonremovable parts. The former is that component of the couplings due to the interaction between the states included in the calculation being performed and that can be removed by appropriate rotation of the adiabatic states into a quasi-diabatic representation. The latter component is because the (infinite number of) other states are being excluded, and this coupling is residual in whichever diabatisation scheme is used. It is these nonremovable couplings that mean eq 19 cannot hold for all pairs of included states and that eq 22 is not an equality. The nonremovable couplings can be safely ignored if the included states form a Hilbert subspace, which is the case if the extended curl and divergence conditions, extensively discussed in ref 12, are (approximately) fulfilled by those states included in the ∂ calculations. Defining F(ijp) = ⟨ψi| ∂p ψj⟩ for some nuclear coordinate p, then, respectively, the extended curl and divergence conditions are ∂F(ijp) ∂q

−

∂F(ijq) ∂p

− [F(ijp), F(ijq)] = 0

∇·Fij = ⟨ψi|∇2 ψj⟩ −

∑ FikFkj k

(26a)

(26b)

where the square brackets denote commutation. These relations hold if the states have negligible coupling to the infinite number of the states in the complement to the subspace.12 A diabatisation scheme that is designed to take the non-removable couplings into account is detailed in ref 55. However, the motivation behind the regularisation and propagation diabatisation methods is to deal with the nonadiabatic elements of quantum nuclear dynamics in a pragmatic manner. As such, the main area of concern is to treat the PES around conical intersections so that wavepacket population transfer can occur in smooth and semirigorous style. As pointed out above, the vector coupling term is singular at such an intersection, that is, it has infinite magnitude there. This infinity, however, is purely due to the interaction of the states involved in the intersection, states that are included in the calculation. The infinite coupling is thus removable and can be dealt with by rotation to the quasi-diabatic basis. As the other, excluded states differ in energy from those at the intersection, the nonremovable component of the couplings is finite and does not prevent one from generating smooth PESs through the intersection. Both the regularisation and propagation diabatisation methods are ways of dealing with the local difficulty in the PES, and as such neither can be considered as truly global diabatisation schemes even though a global phase for the transformation matrix is chosen by the initial conditions of the calculations. However, the propagation diabatisation scheme, unlike the regularisation scheme, reacts to the form of the potential and can remove the necessary couplings if a new intersection is discovered, anywhere in configuration space, during the dynamics calculation. The F


Article


problems discussed above for defining diabatic states, the system is divided into three parts defined by the energy gap between the adiabatic states. In region 1, the adiabatic states are far apart and the removable couplings are small, so to prevent problems due to nonremovable couplings, C is extrapolated from neighboring values. This region is reached when the energy gap between the states exceeds a certain value, set to 2 eV in this study. In this region, the system thus propagates on adiabatic states effectively rotated by a constant matrix. In region 2, the adiabatic states in question are close together, so the removable couplings must now be dealt with by propagating C. This occurs when the state energy gap drops below the threshold for region 1. Considering eq 25, we know C(R) and wish to find C(R + ΔR). To do so we must integrate F over the path R → R + ΔR. The path chosen is simply the straight line in molecular configuration space between the start and end points. We also know F (R) and F (R + ΔR) from the electronic structure calculations. As mentioned in the section “Propagation Diabatisation” we perform a simple, evenly spaced trapezium rule integration along the chosen path. This requires us to know F at each intermediate point, R + kδ, where δ = ΔR/n and 1 ≤ k ≤ (n − 1). To get an accurate integral, we need a fine discretization of the path (we note that a value of n = 40 was found to give good results in the test calculations that follow); however, to calculate F at each point using an electronic structure program would be prohibitively expensive in computational effort. We therefore interpolate F between the start and end points and use those values to perform the necessary integration. As we note from eq 23, each element of F is a ratio of a derivative coupling term and the adiabatic energy gap between the states in question. These two components are both interpolated using Shepard’s method. The adiabatic energies at the point of interest are obtained by using the diabatic energies, gradients, and Hessians at each database point to perform a second-order Taylor expansion, hence getting an approximation of the diabatic energy there. The obtained energies are weighted by the inverse of the square of the Euclidean norm of the difference between the Cartesian coordinates of the molecular geometries of the database point and the point of interest, and then summed. This diabatic energy matrix is then diagonalized to get the adiabatic energies. Subsequently, the adiabatic quantities at the geometry R + ΔR are also weighted and added to provide a contribution to this sum. This gives an unnormalized measure of Vi(R + kδ) (and the same for state j). The derivative coupling in the numerator of eq 23 is also Shepard interpolated with the same weighting scheme, using all database points as well as the new value calculated at R + ΔR. Unlike the energies, where the presence of gradients and Hessians allows a second-order interpolation, these are unavailable for Dij; hence Dij(R + kδ) is just an unnormalized, weighted sum of all available couplings. The two weighted sums can then be divided to give Fij(R + kδ). Note that it is not necessary to normalize either sum because they are both weighted in the same way, and the normalization constants simply cancel. Taking the dot product with δ, we then have the necessary integrand to evaluate the integrals in eq 25. Having generated the integrals, we can then form the matrix exponentials in eq 25 by expanding in a Taylor series. As the matrices are only of low dimensionality (being square of dimension equal to the number of states), and hence their

diabatic energy matrix is simply obtained by the similarity transform in eq 11. The diabatic matrix of gradient vectors, GA, is also found by use of a similarity transform, but this time the sum of the diagonal matrix of adiabatic gradient vectors, GA, and the matrix of the derivative coupling vectors, D (which has zeroes along its diagonal as the couplings are between states only) is transformed.

G̲ D = C†(G̲ A + D̲ )C

(27)

Diabatic Hessian: Powell Update. Neither the calculation of the energies nor of the gradients is problematic as the necessary adiabatic matrix elements are readily available from Gaussian 03 or any other electronic structure program. The problems arise with the calculation of the diabatic Hessian needed for the accurate extrapolation of the PES between database points and for the evaluation of Hamiltonian matrix elements in the nuclear dynamics. However, attempts to transform the adiabatic Hessian to the diabatic representation (both are effectively matrices of matrices, the former being diagonal in terms of the states) using just the transformation matrix, C, revealed that one needs matrix elements of operators that are not available in electronic structure programmes, as far as we are aware. To circumvent this problem, it was decided to use a Hessian update method to get approximate diabatic Hessians at each database point. The method was explained in the recent review66,79 in the case of calculations on a single PES. This used the Powell update formula, which can easily be extended to the multistate case Hij(R + ΔR)=Hij(R) + −

1 (ϵij ⊗ ΔR + ΔR ⊗ ϵij) ΔR·ΔR

ϵij ·ΔR Hij(R) ·ΔR ⊗ ΔR·Hij(R) (ΔR·ΔR)2

(28)

where Hij is the Hessian matrix at the position ij in the diabatic Hessian matrix, HD, and ϵ ij = (G̲ D(R + ΔR) − G̲ D(R))ij is an element (a vector) of the diabatic gradient difference matrix. As the diabatic gradient matrix is already calculated using eq 27, this method is easily implemented and gives the Hessian matrix on each diabatic PES (when i = j) as well as the second-order coupling term (when i ≠ j). From eq 28, it is apparent that we need a reference Hessian to provide a starting point for the update procedure, in our case the diabatic Hessian at the center of the initial wavepacket. As a representation of this it was decided to use the adiabatic Hessian at this point. Even though the adiabatic and diabatic PES are considered equal at this point, the adiabatic and diabatic Hessians will not necessarily be the same. However, we must assume, for this initial study, that they are sufficiently similar that no significant inaccuracy will result in the dynamics. As a side point, the use of the adiabatic Hessian at the reference point means that, for i ≠ j, Hij = 0 there. This means that the commonly used BFGS update formula79 cannot be used as this involves division by products of the Hessian, leading to a singularity in the update procedure (BFGS requires a positive definite initial Hessian). Propagation of the Adiabatic−Diabatic Transformation Matrix. Having considered the initial conditions for the diabatisation scheme, as well as the method for obtaining the diabatic matrix elements from the resulting transformation matrix, we now describe the implementation of the actual propagation of the transformation matrix C. Reflecting the G


Article


phase12,52 and can be seen to arise if we perform the integration in eq 25 along paths around the intersection. The two angles are given by taking clockwise or anticlockwise routes. One solution is simply to choose the added phase to be π/4. This can be implemented in the program by examining each step of the integration in eq 25. By forming the first-order Taylor series of the adiabatic energies of each state around the initial point of the step, it is possible to approximate the adiabatic energies at the end point. If, by using such a projection, it is predicted that the adiabatic energies will change order during the step, that is, the lower-energy state becomes the higher and vice versa, then the phase of the initial transformation matrix simply has π/4 added to give the new matrix so

multiplication is of very low computational cost, this was done to fifth order. Similarly the inversion of the matrix on the lefthand side of that equation is of minimal burden and is achieved using standard LAPACK routines.80 The final point to consider in the solution of eq 25 is the choice of ΔR, the step between database points along which the propagation of the transformation matrix is done. The step is determined by the nuclear dynamics, the GWPs move on the PES, and by comparison to the previously computed database points, we can measure the difference between the new and saved geometries. Once the minimum distance (in terms of the norm of the difference in their Cartesian coordinates) exceeds a certain user-defined threshold (termed dbmin in the program66), a call to the electronic structure program is made to initiate the calculation of a new database point. The difference vector between the new geometry and the closest database geometry is thus ΔR, and the database geometry is R. The only small difference from this general procedure is when the dynamics calculation starts. At this point the GWPs are spread according to the parameters in place for initializing the wave function. For example, if the dynamics is performed using mass-frequency scaled normal modes as a coordinate system, then the GWPs are initially placed at certain, fixed steps along these modes (and along the diagonals between them). As the step lengths are in terms of the normal modes, the difference in the Cartesian coordinates of each position may be greater than twice the threshold distance for calculating a new database point. If this is so, new coordinates are chosen, evenly, between R and R + δ, so that extra electronic structure calculations can be performed and the diabatisation done with a consistently dense set of electronic structure points. In region 3 the adiabatic surfaces cross, and the singularity in the vector coupling term must be dealt with. One of the main aims of any diabatisation scheme is to give PESs that change smoothly through points where the adiabatic PESs are degenerate (and have discontinuous derivatives), for example, at conical intersections. Considering eqs 23 and 24, we see that the integration needed to propagate the transformation matrix runs into problems if the path ΔR passes directly through such a point of degeneracy because the nonadiabatic coupling vector is singular there. Clearly this is a potential problem for our scheme, as the integral is undefined through such points, and indeed was found to be problematic when performing test calculations. In our initial calculations, using only two states, the problematic GWPs were those moving along a point group symmetry plane of the PES to the side of the intersection on which the GWP started. As the transformation matrix C represents a rotation of the adiabatic states by an angle θ ⎛ cos θ −sin θ ⎞ C=⎜ ⎟ ⎝ sin θ cos θ ⎠

⎛1 0⎞ ⎛ 0 −1⎞ ⎟ → C(R + ΔR) = ⎜ ⎟ C(R) = ⎜ ⎝0 1⎠ ⎝1 0 ⎠

(30)

However, on performing test calculations it was found that the diabatic coupling terms that result from such an additional phase are not as smooth as would be expected. A better solution is thus to locate the intersection by projection of the adiabatic energies as described above eq 30, but then to simply extrapolate the diabatic energies and other terms from the database to the end point of the step. The diabatic energy matrix can then be diagonalized, and the generated transformation matrix is then the adiabatic−diabatic transformation matrix that we require as can be seen by examination of eq 11. This point can then be saved to the database, and the transformation matrix can be used for subsequent propagation. In this way smooth couplings are maintained in the vicinity of the conical intersection. Symmetry in the Propagation of the Transformation Matrix. The propagation scheme for the adiabatic−diabatic transformation matrix, described in the above sections, is quite general, but there are a few further subtleties to consider. The first is the symmetry of the PES generated. Consider Figure 1, which shows the initial positions of 25 GWPs along two massfrequency scaled normal modes of butatriene. The origin of the coordinate system is the Franck−Condon point of the molecule, which has D2h symmetry. The x-axis is the stretch of the central carbon−carbon bond, labeled 14Ag, and the y-axis is the torsion-type mode of the linear carbon skeleton, labeled 5Au. Clearly, by symmetry considerations, the PES should be symmetric along the 5Au mode (i.e., it has a plane of symmetry along the x-axis). As described in the previous sections, the resultant PES would not maintain this symmetry, even with the symmetrization of the wavepacket outlined above. The reason for the breaking of symmetry can be explained as follows: GWP 1 is at the reference geometry, the center of the wavepacket, and the PES at the center of GWP 2 gains a diabatic representation by solution of eq 25 starting at the origin. As the database has but a single entry, corresponding to GWP 1, the propagated F is extrapolated from this point only. The transformation matrix at the center of GWP 3 is likewise gained by propagation from GWP 1, but this time it uses information from two database points to get F along the path. The same is repeated for GWP 4, but with F interpolated from the three previous database points. The transformation matrix at GWP 4 will thus be different than that at GWP 2; however, GWPs 2 and 4 are centered at symmetry equivalent geometries, so the matrices

(29)

we can take θ = 0 when the transformation matrix is the identity matrix. The problem we have is to define θ to the far side of the intersection, where the diabatic states have crossed (on the initial side, adiabatic state 1 corresponds to diabatic state 1, but it corresponds to diabatic state 2 on the other). Considering eqs 11 and 27, it is clear that a choice of θ = π/4 or −π/4 leads to the same diabatic energies and gradients (and by extension the Hessian that is constructed from the gradients) because C(π/4) = −C(−π/4), and the sign cancels in a similarity transform. This uncertainty relates to the Berry H


Article

The Journal of Physical Chemistry A should be the adjoint of one another. This will not be the case; hence, the diabatic surfaces will lose symmetry. To overcome this problem we consider the GWPs in terms of the shells to which they belong. Only database points corresponding to GWPs in shells closer to the reference geometry than they are, are used in the generation of F, and as a result, each GWP in a particular shell is treated on an equal footing. In the current example, only the database point corresponding to GWP 1 is used to diabatise the PES at GWPs 2, 3, 4, and 5. Likewise, from the diagram, the PES centered at GWPs 6 to 9 inclusive will be diabatised using information from database points 1 to 5. The same procedure is maintained during the dynamics calculation to ensure that symmetry is maintained. Only database points generated at previous time steps and those from inner shell GWPs at the same time step are used. A similar concern arises for the Powell updated Hessian as described in the section Diabatic Hessian: Powell Update, and the problem is solved in a similar way. At the end of each set of calculations of new database points and transformation matrices, the Hessians for all database points are regenerated using the very first as reference. Only database points strictly closer to the reference geometry than the point under consideration are used in generating the Hessian. This ensures that symmetry equivalent points do not contribute to the Hessians of one another, and are generated using the same set of data.

Figure 2. Pictorial representation of the normal modes of butatriene used as coordinates in this work. (a) 5Au: molecular torsion. (b) 8Ag: symmetric outer C−C bond stretching. (c) 12Ag: symmetric CH2 bending. (d) 14Ag: symmetric stretching of the central C−C bond. (e) 15Ag: symmetric C−H bond stretching.

regularized diabatisation calculations) at the GWP centers to build the database needed to obtain the diabatic PES. New database points were calculated once a GWP had moved so that it was more than 0.1 bohr from the nearest database point (measured by taking the norm of the Cartesian vector, ΔR). The dynamics was performed using a fifth-order Runge−Kutta integrator with accuracy cutoff of 1 × 10−6 to propagate both the A-coefficients and the GWP parameters. The propagations were allowed to run for 100 fs with data output every 0.5 fs. In Figure 3a we plot the adiabatic PES generated during the course of the calculation using the propagation diabatisation scheme. In total 153 database points were generated during the 100 fs of the propagation. It is clear from this plot that a conical intersection is present at ∼1.7 units along the 14Ag mode and at 0 along the 5Au (i.e., at a planar geometry). In Figure 3b we present the diabatic surfaces resulting from the propagation diabatisation (using the set of points from Figure 3a). It is apparent from this that we have generated a pair of smoothly varying, quasi-diabatic surfaces with a crossing coinciding with the observed conical intersection. This can be compared to Figure 3c, which shows the quasi-diabatic states generated using the regularisation diabatisation scheme (from 154 adiabatic points in this case). The similarities between Figure 3b,c are clear; the shapes of the surfaces are essentially the same with the crucial crossing being at the same point. Recall that the regularisation diabatisation relies on locating the conical intersection with a prior electronic structure calculation (which optimizes to a geometry exhibiting a degeneracy), so the crossing point was determined prior to the dynamics being performed. There is no such initial information needed with the propagation diabatisation scheme; the crossing is found automatically. That the crossings match between the two different calculations is validation of the method that we have proposed in this work. In comparison with the adiabatic surfaces in Figure 3a we see that the propagation diabatisation diabatic surfaces better represent the higher-energy parts at large negative values of the 14Ag coordinate than do the surfaces generated by the regularisation diabatisation method, indicating the greater degree of flexibility in the new scheme. This improvement is a reflection of the division of the configuration space in the propagation diabatisation scheme into regions as defined in the section Propagation of the Adiabatic−Diabatic Transformation Matrix, these portions of the PES being separated by more than

■

RESULTS Two-Mode Butatriene. To test the propagation diabatisation scheme we performed a series of calculations on the singly charged cation of butatriene. Butatriene is a classic system dominated by vibronic coupling and was originally used to validate the implementation of the regularisation diabatisation scheme into DD-vMCG.65 The initial DD-vMCG calculations were performed using a two-dimensional subset of the 18 mass-frequency scaled normal modes that describe the vibrational motions of butatriene. These were the 5Au and 14Ag modes obtained from a CASSCF(6,6)/3-21G* calculation, using Gaussian 03, to optimize the neutral ground-state geometry, followed by a diagonalization of the molecular Hessian to get the modes. The frequencies of the modes are 767.63 and 2196.18 cm−1, respectively, representing the torsion of the carbon backbone and the symmetric stretching of the central carbon−carbon double bond (as seen in Figure 2). In addition a state-averaged CASSCF(5,6)/3-21G* calculation was performed on the cation ground and first-excited states (weighted equally) to locate the position of the conical intersection needed as a reference for the regularized diabatisation calculations. DD-vMCG nuclear dynamics calculations were then performed on butatriene, once with the regularized diabatisation scheme and once with the new, propagation diabatisation method. Dynamics proceeded on the ground and first-excited states of the cation only. The initial wavepacket was a Gaussian function of width 1/√2 and zero momentum along both modes, centered at the Franck−Condon point (the origin of the normal mode coordinates) and placed on the excited-state PES. The basis set used consisted of a single set of 25 GWPs of frozen width, 1/√2, initially centered as illustrated in Figure 1 (the spacings are multiples of 0.649 units). The state-averaged CASSCF(5,6)/3-21G* method was used to generate the adiabatic energies and gradients (and Hessians for the I


Article


Figure 3. PESs for the butatriene cation along the 5Au and 14Ag mass-frequency scaled normal modes. Energy and coordinate origins represent the minimum energy structure in the neutral ground state. (a) Adiabatic surfaces, D0 and D1, calculated at the state-averaged CASSCF(5,6)/3-21G* level using Gaussian 03. (b) Quasi-diabatic surfaces generated from those in (a) using the new propagation diabatisation scheme outlined in the section Propagation Diabatisation. (c) Quasi-diabatic surfaces generated using the regularisation diabatisation scheme outlined in the section Regularisation Diabatisation. A different database of adiabatic points than those in (a) and (b) was used due to the different dynamics in this case. (d) Diabatic coupling terms generated by the propagation diabatisation scheme (i.e., the matrix element VD12 defined in eqs 11 and 9) as a function of molecular geometry.

more complete description of the butatriene cation would require consideration of at least three states. Having seen that the new diabatisation scheme produces PES in good agreement with the established method, we need to see whether the actual nuclear dynamics correspond or not. To do so, we present a plot in Figure 4 of the wavepacket population

2 eV and hence in region 1, the pseudoadiabatic part. The plot of the diabatic coupling terms in Figure 3d reflects this division. At positive values of the 14Ag coordinate the couplings vary linearly with changes in the 5Au coordinate. The slope along the nontotally symmetric mode varies as a function of the 14Ag mode. At negative values of the totally symmetric mode the linear relationship is maintained around the planar molecular geometry (5Au = 0), but there is a distinct deviation at larger absolute values of the coordinate, corresponding to region 1. Here the transformation matrix is no longer being propagated, so the couplings deviate from their previous pattern. There is also a small difference between the regularisation and propagation diabatisation generated surfaces to the right of the conical intersection (where 14Ag > 2) at geometries close to planar. This difference is a reflection of the different methods of diabatisation around the intersection, the regularisation diabatisation surfaces are expanded around that point, while the propagation diabatisation surfaces are generated by GWPs that have passed through region 3, and particularly directly through the intersection. The surfaces further away from planar, and hence from the intersection, are more similar between the two methods, any differences having been smoothed out. Another feature that was noticed during the test calculations was that the vector coupling terms suddenly became orthogonal to the motion of the GWPs at certain geometries far from the intersection, whereas, closer in, they are well-aligned. Further investigation revealed the cause of this to be the presence of a third, intruder state coming down in energy and crossing the upper state in the model. The horseshoe-shaped region, where this is the case, is at energies too high to have a significant effect on the dynamics, however, as it was sampled by GWPs that moved to that region because of the spreading of the wavepacket on the ground state. We can thus conclude that a

Figure 4. Wavepacket population of the first-excited, diabatic, cation state of butatriene as a function of propagation time using a twodimensional model consisting of the 5Au and 14Ag normal modes. Solid line corresponds to propagation using the MCTDH method on surfaces fitted, using the vibronic coupling Hamiltonian, to electronic energies calculated at the CASSCF(5,6)/3-21G* level. Dashed line corresponds to DD-vMCG calculations using the regularized diabatisation method. Dotted line corresponds to DD-vMCG calculations using the new, propagation diabatisation method. Both DD-vMCG calculations also used state-averaged CASSCF(5,6)/321G* electronic energies. J


Article

The Journal of Physical Chemistry A of the initially populated, diabatic state (the excited state at the Franck−Condon point) for DD-vMCG using both diabatisation schemes as well as for an exact wavepacket dynamics calculation on PESs fitted to adiabatic points calculated along the two modes in question28 using the vibronic coupling Hamiltonian approach. This calculation was also performed using the Quantics package. The similarity between the three plots is clear. Starting from a fully populated upper state, there is rapid depopulation over the first 15 fs or so (down to the 20−40% range) as the wavepacket encounters the surface crossing seam and moves onto the other diabatic surface. This depopulation is followed by a repopulation to ∼50−60% over the next 10 fs with a subsequent depopulation to 15−25%. Beyond 40 fs the two DD-vMCG calculations give rapidly oscillating populations around a midpoint of ∼40%. The full calculation on the fitted surface shows clearer structure for these later times with more obvious periods of increasing and decreasing populations, but overall around an average of ∼50% of the initial population. It should be pointed out that we do not expect the DDvMCG plots to converge on the exact dynamics plot, even given an infinite basis of GWPs and more densely spaced database points, because the PESs we are comparing are constructed differently. For the exact dynamics, the surfaces are a best fit to a given analytic form expanded around the Franck− Condon point. For DD-vMCG, the nuclear wavepacket moves on PESs that are interpolated between the given calculated points, and as such the PESs match the calculated energies exactly at the database geometries. For a fitted surface, this is only generally true at the central point of expansion. For this reason and because there is no initial assumption made about the functional form of the surfaces then, in principle, the DDvMCG calculations should give a more accurate representation of the PESs. It is, however, reassuring and a key test that dynamics calculations using two different types of method give results in very good qualitative agreement. The main comparison we need to make to validate the new diabatisation method is between the two DD-vMCG plots in Figure 4, which use the different diabatisation schemes. In this respect we see that the propagation diabatisation scheme is working well, the qualitative shape of the plot matches well to that generated by the dynamics using the older method. There are quantitative differences, however: in particular the level of depopulation is greater for the regularisation diabatisation calculation during the first 40 fs, although the difference between the two plots decreases over this time. The initial decrease in population is to ∼21% of the initial value for the regularisation diabatisation method and to just under 37% for the propagation diabatisation calculation, while the first recovery is to ∼52% and 62%, respectively. The decreasing difference is followed through to the final major depopulation where the state populations fall to 13% and 17% for the regularisation and propagation diabatisation methods. Beyond 40 fs the plots continue to follow one another qualitatively despite the more rapid oscillations; for example, the broader maximum and minimum between 80 and 90 fs is seen in both plots. We thus see that the propagation diabatisation method is more than capable of successfully reproducing the results of the established method in this case. Five-Mode Butatriene. In Figure 5 we show the equivalent plots of diabatic state populations for the butatriene cation, as seen in the preceding section, but this time using a five-mode model. The extra modes included are those totally symmetric

Figure 5. Wavepacket population of the first-excited, diabatic, cation state of butatriene as a function of propagation time using a fivedimensional model consisting of the 5Au, 8Ag, 12Ag, 14Ag, and 15Ag normal modes. Solid line corresponds to propagation using the MCTDH method on surfaces fitted, using the vibronic coupling Hamiltonian, to electronic energies calculated at the CASSCF(5,6)/321G* level. Dashed line corresponds to DD-vMCG calculations using the regularized diabatisation method. Dotted line corresponds to DDvMCG calculations using the new, propagation diabatisation method. Both DD-vMCG calculations also used state-averaged CASSCF(5,6)/ 3-21G* electronic energies.

ones labeled 8Ag (stretching of the outer C−C bonds), 12Ag (CH2 bending), and 15Ag (CH stretching). All of the modes can be seen in Figure 2. Once again we compare calculations on the fitted surfaces,28 this time using the MCTDH algorithm81 rather than propagation on the full grid. The multiset formulation (where different basis functions are used for each electronic state) was used with eight two-dimensional, timedependent basis functions per state for the pairs of modes 5Ag/ 14Ag and 12Ag/15Ag and six one-dimensional functions for each state along the 8Ag mode. The DD-vMCG calculations were performed using 51 five-dimensional GWPs in the single set formulation, this number being used to ensure a symmetric spread of the basis functions along each mode. From this figure we can see that the DD-vMCG populations follow those generated by the MCTDH calculation quite well. We get a similar picture of rapid depopulation of the upper state followed by partial repopulation. Beyond ∼30 fs the two methods no longer follow each other, reflecting the different methods used to construct the PESs. However, once again we see good agreement between the regularisation and propagation diabatisation results indicating the utility of the propagation diabatisation method. One further point that should be remarked upon is the relatively jagged change of the population that is seen with the DD-vMCG calculations, both in the two and five-dimensional situations, especially when compared to the smooth shape of the plots generated by the calculations on the fitted surfaces. We generally see this type of behavior when performing vMCG calculations and attribute it to incompleteness of the GWP basis set and also a certain irregularity in the motion of those GWPs, the latter in part caused by the symmetrization procedure outlined in the section DD-vMCG. Work is ongoing in our group to alleviate these more general issues with propagation scheme. K


Article


■

(4) Köppel, H.; Cederbaum, L. S.; Domcke, W. Strong non-adiabatic effects and conical intersections in molecular-spectroscopy and unimolecular decay - C2H4+. J. Chem. Phys. 1982, 77, 2014−2022. (5) Köppel, H.; Domcke, W.; Cederbaum, L. S. Multimode molecular Dynamics beyond the Born-Oppenheimer approximation. Advances in Chemical Physics 1984, 57, 59−246. (6) Manthe, U.; Köppel, H. Dynamics on potential energy surfaces with a conical intersection: Adiabatic, intermediate and diabatic behaviour. J. Chem. Phys. 1990, 93, 1658−1669. (7) Lengsfield, B. H.; Yarkony, D. R. Non-adiabatic interactions between potential energy surfaces: Theory and applications. Advances in Chemical Physics 1992, 82, 1−71. (8) Yarkony, D. R. Diabolical conical intersections. Rev. Mod. Phys. 1996, 68, 985−1013. (9) Yarkony, D. R. Conical intersections: diabolical and often misunderstood. Acc. Chem. Res. 1998, 31, 511−518. (10) Baer, M. Introduction to the theory of electronic non-adiabatic coupling terms in molecular systems. Phys. Rep. 2002, 358, 75−142. (11) Mebel, A. M.; Yahalom, A.; Englman, R.; Baer, M. The study of conical intersections between consecutive pairs of the five lowest 2A′ states of the C2H molecule. J. Chem. Phys. 2001, 115, 3673−3689. (12) Baer, M. Beyond Born-Oppenheimer Electronic Nonadiabatic Coupling Terms and Conical Intersections; John Wiley and Sons Inc: Hoboken, NJ, 2006. (13) Cederbaum, L. S.; Köppel, H.; Domcke, W. Multimode vibronic coupling effects in molecules. Int. J. Quantum Chem. 1981, 20, 251− 267. (14) Neville, S.; Worth, G. A reinterpretation of the electronic spectrum of pyrrole: A quantum dynamics study. J. Chem. Phys. 2014, 140, 034317. (15) Giri, K.; Chapman, E.; Sanz, C. S.; Worth, G. A full-dimensional coupled-surface study of the photodissociation dynamics of ammonia using the multiconfiguration time-dependent Hartree method. J. Chem. Phys. 2011, 135, 044311. (16) Li, Z. H.; Valero, R.; Truhlar, D. G. Improved direct diabatization and coupled potential energy surfaces for the photodissociation of ammonia. Theor. Chem. Acc. 2007, 118, 9−24. (17) Worth, G. A.; Robb, M. A.; Lasorne, B. Solving the timedependent Schrödinger Equation in one step: Direct Dynamics of Non-adiabatic Systems. Mol. Phys. 2008, 106, 2077−2091. (18) Persico, M.; Granucci, G. An overview of nonadiabatic dynamics simulation methods, with focus on the direct approach versus the fitting of potential energy surfaces. Theor. Chem. Acc. 2014, 133, 1526. (19) Tully, J. C.; Preston, R. K. Trajectory surface hopping approach to nonadiabatic molecular collisions: The reaction of H+ with D2. J. Chem. Phys. 1971, 55, 562−572. (20) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061−1071. (21) Tully, J. C. Mixed quantum-classical dynamics. Faraday Discuss. 1998, 110, 407−419. (22) Coker, D. F. Computer Simulation Methods for Nonadiabatic Dynamics in Condensed Systems. Computer Simulation in Chemical Physics. Dordrecht 1993, 315−377. (23) Coker, D. F.; Xiao, L. Methods for molecular dynamics with nonadiabatic transitions. J. Chem. Phys. 1995, 102, 496−510. (24) Topaler, M. S.; Hack, M. D.; Allison, T. C.; Liu, Y.-P.; Mielke, S. L.; Schwenke, D. W.; Truhlar, D. G. Validation of trajectory surface hopping methods against accurate quantum mechanical dynamics and semiclassical analysis of electronic-to-vibrational energy transfer. J. Chem. Phys. 1997, 106, 8699−8709. (25) Hack, M. D.; Truhlar, D. G. Nonadiabatic Trajectories at an Exhibition. J. Phys. Chem. A 2000, 104, 7917−7926. (26) Richter, M.; Marquetand, P.; González-Vázquez, J.; Sola, I.; González, L. SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. J. Chem. Theory Comput. 2011, 7, 1253−1258. (27) Worth, G. A.; Robb, M. A. Applying direct molecular dynamics to non-adiabatic systems. Adv. Chem. Phys. 2002, 124, 355−432.

CONCLUSIONS We have presented an implementation of a new diabatisation scheme for the direct dynamics version of the vMCG algorithm for fully quantum nuclear dynamics. The method is based upon the propagation of the adiabatic−diabatic transformation matrix by use of the nonadiabatic vector coupling terms obtainable from electronic structure packages. The implementation has many technical subtleties, as outlined in the section so titled, but these have resulted in a method that should be generally applicable to any molecular system as well as relatively simple to use for the nonspecialist user. Simplicity and generality are the key benefits over other, more involved schemes. In our initial, test calculations, the method has been shown to give results that compare well to the already established regularisation diabatisation scheme. Our method is in the same spirit as regularisation diabatisation, diabatising the PESs around conical intersections by dealing with the removable parts of the nonadiabatic coupling terms in a smooth fashion. The propagation diabatisation method overcomes some of the weaknesses of the older scheme though. In particular, there is no need to locate any conical intersections prior to performing the dynamics. As the method is not defined by the location of a specific intersection, it should work equally well whenever an intersection is located during the course of the dynamics. It will be the focus of upcoming work to test the propagation diabatisation method on systems exhibiting nontrivial conical intersection seams, as well as those with multiple, unconnected intersections. The propagation diabatisation method can also, in principle, diabatise any number of electronic states as the vector coupling between all pairs of the states can be calculated (i.e., all elements of F can be evaluated). Testing of this feature of the method will also be part of our upcoming work, the first target being the butatriene system used in this work where there is a third, intruder state coming down to cross the first excited state at geometries relatively far from the Franck−Condon point. Inclusion of more than two states will be a key test of the utility of our method.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was funded by the EPSRC as part of the development of the QUANTICS quantum dynamics package EP/ K037943. The program is available on request from the authors. All data presented in this paper is also available on request.

■

REFERENCES

(1) Worth, G. A.; Cederbaum, L. S. Beyond Born-Oppenheimer: Conical intersections and their impact on molecular dynamics. Annu. Rev. Phys. Chem. 2004, 55, 127−158. (2) Conical intersections: Electronic structure, dynamics and spectroscopy; Domcke, W., Yarkony, D. R., Köppel, H., Eds.; World Scientific: Singapore, 2004. (3) Domcke, W.; Köppel, H.; Cederbaum, L. S. Spectroscopic effects of conical intersections of molecular potential surfaces. Mol. Phys. 1981, 43, 851−875. L


Article


using frozen-width variational Gaussian product basis functions. J. Chem. Phys. 2012, 137, 22A548. (48) Allan, C.; Lasorne, B.; Worth, G.; Robb, M. A Straightforward Method of Analysis for Direct Quantum Dynamics: Application to the Photochemistry of a Model Cyanine. J. Phys. Chem. A 2010, 114, 8713−8729. (49) Asturiol, D.; Lasorne, B.; Worth, G.; Blancafort, L.; Robb, M. Exploring the sloped-to-peaked S2/S1 seam of intersection of thymine with electronic structure and direct quantum dynamics calculations. Phys. Chem. Chem. Phys. 2010, 12, 4949−4958. (50) Araújo, M.; Lasorne, B.; Magalhaes, A.; Bearpark, M.; Robb, M. Controlling Product Selection in the Photodissociation of Formaldehyde: Direct Quantum Dynamics from the S1 Barrier. J. Phys. Chem. A 2010, 114, 12016−12020. (51) Araújo, M.; Lasorne, B.; Magalhaes, A.; Worth, G.; Bearpark, M.; Robb, M. The molecular dissociation of formaldehyde at medium photoexcitation energies: A quantum chemistry and direct quantum dynamics study. J. Chem. Phys. 2009, 131, 144301. (52) Tannor, D. J. Introduction to Quantum Mechanics A Time Dependent Perspective; University Science Books: South Orange, NJ, 2007. (53) Baer, M. Adiabatic and diabatic representations for atommolecule collisions: Treatment of the collinear arrangement. Chem. Phys. Lett. 1975, 35, 112−118. (54) Mead, C. A.; Truhlar, D. G. Conditions for the definition of a strictly diabatic electronic basis for molecular systems. J. Chem. Phys. 1982, 77, 6090−6098. (55) Evenhuis, C.; Martínez, T. A scheme to interpolate potential energy surfaces and derivative coupling vectors without performing a global diabatization. J. Chem. Phys. 2011, 135, 224110. (56) Godsi, O.; Evenhuis, C.; Collins, M. Interpolation of multidimensional diabatic potential energy matrices. J. Chem. Phys. 2006, 125, 104105. (57) Evenhuis, C.; Lin, X.; Zhang, D.; Yarkony, D.; Collins, M. interpolation of diabatic potential-energy surfaces: Quantum dynamics on ab initio surfaces. J. Chem. Phys. 2005, 123, 134110. (58) Alijah, A.; Baer, M. The Electronic Adiabatic-Diabatic Transformation Matrix: A Theoretical and Numerical Study of a Three-State System. J. Phys. Chem. A 2000, 104, 389−396. (59) Cave, R.; Stanton, J. Block diagonalization of the equation-ofmotion coupled cluster effective Hamiltonian: Treatment of diabatic potential constants and triple excitations. J. Chem. Phys. 2014, 140, 214112. (60) Hoyer, C.; Xu, X.; Ma, D.; Gagliardi, L.; Truhlar, D. Diabatization based on the dipole and quadrupole: The DQ method. J. Chem. Phys. 2014, 141, 114104. (61) Nakamura, H.; Truhlar, D. The direct calculation of diabatic states based on configurational uniformity. J. Chem. Phys. 2001, 115, 10353−10372. (62) Nakamura, H.; Truhlar, D. Direct diabatization of electronic states by the fourfold way. II. Dynamical correlation and rearrangement processes. J. Chem. Phys. 2002, 117, 5576−5593. (63) Nakamura, H.; Truhlar, D. Extension of the fourfold way for calculation of global diabatic potential energy surfaces of complex, multiarrangement, non Born-Oppenheimer systems: Aplication to HNCO(S0,S1). J. Chem. Phys. 2003, 118, 6816−6829. (64) Worth, G. A.; Meyer, H.-D.; Beck, M. H.; Burghardt, I.; Jäckle, A.; Burghardt, I.; Giri, K.; Lasorne, B.; Otto, F.; Richings, G. W.; Schröder, M. Quantics: A Computer Package for Quantum Dynamics, Development Version 1.0. 2015; See http://stchem.bham.ac.uk/ quantics. (65) Worth, G. A.; Robb, M. A.; Burghardt, I. Non-adiabatic Direct Dynamics using Variational Gaussian Wavepackets: The X̃ /Ã Manifold of the Butatriene Cation. Faraday Discuss. 2004, 127, 307−323. (66) Richings, G.; Polyak, I.; Spinlove, K.; Worth, G.; Burghardt, I.; Lasorne, B. Quantum Dynamics Simulations using Gaussian Wavepacket: The vMCG Method. Int. Rev. Phys. Chem. 2015, 34, 269−308. (67) Ischtwan, J.; Collins, M. A. Molecular potential energy surfaces by interpolation. J. Chem. Phys. 1994, 100, 8080−8088.

(28) Worth, G.; Hunt, P.; Robb, M. Non-adiabatic dynamics: A comparison of surface hopping direct dynamics with quantum wavepacket calculations. J. Phys. Chem. A 2003, 107, 621−631. (29) Smith, B. R.; Bearpark, M. J.; Robb, M. A.; Bernardi, F.; Olivucci, M. ‘Classical wavepacket’ dynamics through a conical intersection. Application to the S1/S0 photochemistry of benzene. Chem. Phys. Lett. 1995, 242, 27−32. (30) Bearpark, M. J.; Bernardi, F.; Olivucci, M.; Robb, M. A.; Smith, B. R. Can fulvene S1 decay be controlled? A CASSCF study with MMVB dynamics. J. Am. Chem. Soc. 1996, 118, 5254−5260. (31) Bearpark, M. J.; Bernardi, F.; Clifford, S.; Olivucci, M.; Robb, M. A.; Smith, B. R.; Vreven, T. The azulene S1 state decays via a conical intersection: A CASSCF study with MMVB dynamics. J. Am. Chem. Soc. 1996, 118, 169−175. (32) Richter, M.; Marquetand, P.; González-Vázquez, J.; Sola, I.; González, L. Femtosecond Intersystem Crossing in the DNA Nucleobase Cytosine. J. Phys. Chem. Lett. 2012, 3, 3090−3095. (33) Mitrić, R.; Petersen, J.; Bonacić-Koutecký, V. Laser-field-induced surface-hopping method for the simulation and control of ultrafast photodynamics. Phys. Rev. A: At., Mol., Opt. Phys. 2009, 79, 053416. (34) Mitrić, R.; Petersen, J.; Wohlgemuth, M.; Werner, U.; BonacićKoutecký, V. Field-induced surface hopping method for probing transition state nonadiabatic dynamics of Ag3. Phys. Chem. Chem. Phys. 2011, 13, 8690−8696. (35) Mitrić, R.; Petersen, J.; Wohlgemuth, M.; Werner, U.; BonacićKoutecký, V.; Wöste, L.; Jortner, J. Time-Resolved Femtosecond Photoelectron Spectroscopy by Field-Induced Surface Hopping. J. Phys. Chem. A 2011, 115, 3755−3765. (36) Lisinetskaya, P.; Mitrić, R. Simulation of laser-induced coupled electron-nuclear dynamics and time-resolved harmonic spectra in complex systems. Phys. Rev. A: At., Mol., Opt. Phys. 2011, 83, 033408. (37) Ben-Nun, M.; Martínez, T. J. Ab Initio Quantum Molecular Dynamics. Adv. Chem. Phys. 2002, 121, 439−512. (38) Hudock, H. R.; Levine, B. G.; Thompson, A. L.; Satzger, H.; Townsend, D.; Gador, N.; Ullrich, S.; Stolow, A.; Martínez, T. J. Ab initio molecular dynamics and time-resolved photoelectron spectroscopy of electronically excited uracil and thymine. J. Phys. Chem. A 2007, 111, 8500−8508. (39) Mori, T.; Glover, W.; Schuurman, M.; Martínez, T. Role of Rydberg States in the Photochemical Dynamics of Ethylene. J. Phys. Chem. A 2012, 116, 2808−2818. (40) Ko, C.; Levine, B.; Toniolo, A.; Manohar, L.; Olsen, S.; Werner, H.-J.; Martínez, T. J. Ab initio excited-state dynamics of the photoactive yellow protein chromophore. J. Am. Chem. Soc. 2003, 125, 12710−12711. (41) Martínez, T. J.; Ben-Nun, M.; Levine, R. D. Multi-electronicstate molecular dynamics: A wave function approach with applications. J. Phys. Chem. 1996, 100, 7884−7895. (42) Ben-Nun, M.; Martínez, T. J. Nonadiabatic molecular dynamics: Validation of the multiple spawning method for a multidimensional problem. J. Chem. Phys. 1998, 108, 7244−7257. (43) Lasorne, B.; Sicilia, F.; Bearpark, M. J.; Robb, M. A.; Worth, G. A.; Blancafort, L. Automatic generation of active coordinates for quantum dynamics calculations: Application to the dynamics of benzene photochemistry. J. Chem. Phys. 2008, 128, 124307−10. (44) Lasorne, B.; Bearpark, M. J.; Robb, M. A.; Worth, G. A. Controlling S1/S0 Decay and the Balance between Photochemistry and Photostability in Benzene: A Direct Quantum Dynamics Study. J. Phys. Chem. A 2008, 112, 13017−13027. (45) Lasorne, B.; Robb, M. A.; Worth, G. A. Direct quantum dynamics using variational multi-configuration Gaussian wavepackets. Implementation details and test case. Phys. Chem. Chem. Phys. 2007, 9, 3210−3227. (46) Mendive-Tapia, D.; Lasorne, B.; Worth, G.; Bearpark, M.; Robb, M. Controlling the mechanism of fulvene S1/S0decay: switching off the stepwise population transfer. Phys. Chem. Chem. Phys. 2010, 12, 15725−15733. (47) Mendive-Tapia, D.; Lasorne, B.; Worth, G.; Robb, M.; Bearpark, M. Towards converging non-adiabatic direct dynamics calculations M


Article

The Journal of Physical Chemistry A (68) Collins, M. A.; Zhang, D. H. Application of interpolated potential energy surfaces to quantum reactive scattering. J. Chem. Phys. 1999, 111, 9924−9931. (69) Frankcombe, T.; Collins, M.; Worth, G. Converged quantum dynamics with modified Shepard interpolation and Gaussian wave packets. Chem. Phys. Lett. 2010, 489, 242−247. (70) Thiel, A.; Köppel, H. Proposal and numerical test of a simple diabatization scheme. J. Chem. Phys. 1999, 110, 9371−9383. (71) Köppel, H.; Gronki, J.; Mahapatra, S. Construction scheme for regularized diabatic states. J. Chem. Phys. 2001, 115, 2377−2388. (72) Köppel, H.; Schubert, B. The concept of regularized diabatic states for a general conical intersection. Mol. Phys. 2006, 104, 1069− 1079. (73) Frisch, M. J. et al. Gaussian 03, Revision C.02; Gaussian, Inc: Pittsburgh, PA, 2003. (74) Esry, B.; Sadeghpour, H. Split diabatic representation. Phys. Rev. A: At., Mol., Opt. Phys. 2003, 68, 042706. (75) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. J. Chem. Phys. 2001, 114, 10608−10615. (76) Mai, S.; Marquetand, P.; González, L. Non-adiabatic and intersystem crossing dynamics in S02. II. The role of triplet states in the bound state dynamics studied by surface-hopping simulations. J. Chem. Phys. 2014, 140, 204302. (77) Mai, S.; Marquetand, P.; González, L. A General Method to Describe Intersystem Crossing Dynamics in Trajectory Surface Hopping. Int. J. Quantum Chem. 2015, 115, 1215−1231. (78) Worth, G. A.; Beck, M. H.; Jäckle, A.; Meyer, H.-D. The MCTDH Package, Version 8.3, See http://www.pci.uni-heidelberg.de/ tc/usr/mctdh/. 2003. (79) Frankcombe, T. Using Hessian update formulae to construct modified Shepard interpolated potential energy surfaces: application to vibrating surface atoms. J. Chem. Phys. 2014, 140, 114108. (80) Anderson, E.; Bai, Z.; Bischof, C.; Blackford, S.; Demmel, J.; Dongarra, J.; du Croz, J.; Greenbaum, A.; Hammarling, S.; McKenney, A.; Sorensen, D. LAPACK Users’ Guide, 3rd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1999. (81) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The multiconfiguration time-dependent Hartree method: A highly efficient algorithm for propagating wavepackets. Phys. Rep. 2000, 324, 1−105.

N


A Practical Diabatisation Scheme for Use with the Direct-Dynamics

Recommend Documents