Construction of Diabatic Hamiltonian Matrix from ab Initio Calculated

Article pubs.acs.org/JPCA

Construction of Diabatic Hamiltonian Matrix from ab Initio Calculated Molecular Symmetry Adapted Nonadiabatic Coupling Terms and Nuclear Dynamics for the Excited States of Na3 Cluster Saikat Mukherjee,† Sudip Bandyopadhyay,‡ Amit Kumar Paul,† and Satrajit Adhikari*,† †

Department of Physical Chemistry and Raman Center for Atom, Molecule and Optical Sciences, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700 032, India ‡ Department of Chemistry, Maulana Azad College, Kolkata 700 013, India S Supporting Information *

ABSTRACT: We present the molecular symmetry (MS) adapted treatment of nonadiabatic coupling terms (NACTs) for the excited electronic states (22E′ and 12A1′ ) of Na3 cluster, where the adiabatic potential energy surfaces (PESs) and the NACTs are calculated at the MRCI level by using an ab initio quantum chemistry package (MOLPRO). The signs of the NACTs at each point of the configuration space (CS) are determined by employing appropriate irreducible representations (IREPs) arising due to MS group, and such terms are incorporated into the adiabatic to diabatic transformation (ADT) equations to obtain the ADT angles. Since those sign corrected NACTs and the corresponding ADT angles demonstrate the validity of curl condition for the existence of three-state (22E′ and 12A1′ ) sub-Hilbert space, it becomes possible to construct the continuous, single-valued, symmetric, and smooth 3 × 3 diabatic Hamiltonian matrix. Finally, nuclear dynamics has been carried out on such diabatic surfaces to explore whether our MS-based treatment of diabatization can reproduce the pattern of the experimental spectrum for system B of Na3 cluster.

I. INTRODUCTION The Born−Oppenheimer (BO) treatment1,2 introduces the distinction between the fast and slow moving electrons and nuclei with two important outcomes, viz. the adiabatic potential energy surfaces (PESs) and the nonadiabatic coupling terms (NACTs). When a molecular process is assumed to occur exclusively on the ground adiabatic surface, the relevant Schrödinger equation (SE) is expected to provide enough accurate solutions for the observables, such as reactive/ nonreactive cross-sections or spectroscopic quantities. On the contrary, the excited electronic states affect the ground state dynamics very strongly due to the so-called nonadiabatic coupling terms, and the BO approximate equation fails to calculate correct transition probabilities. Calculations on charge transfer3−7 reactions and scattering cross-sections8−19 on single-, two-, and three-surfaces show serious discrepancies between single and multiple surface results. While surrounding a point of degeneracy in the configuration space (CS), Longuet−Higgins (LH)20,21 and colleagues22,23 observed that the BO adiabatic eigenfunction may acquire a phase and flip their sign. Varandas, Tennyson, and Murrel applied this LH theorem on a realistic system LiNaK to demonstrate24 that the dominant coefficient, ci, of the electronic eigenfunction actually undergoes a sign change along a closed © 2013 American Chemical Society

loop encircling a conical intersection (CI). Though Herzberg and Longuet−Higgins (HLH)25 predicted the multivaluedness of the electronic eigenfunction for the Jahn−Teller (JT) CI model and corrected it in an ad hoc manner, it had an important consequence, namely, such a model Hamiltonian under transformation produces nonadiabatic coupling terms between the two adiabatic potential energy surfaces. On the contrary, the adiabatic-to-diabatic transformation (ADT) technique to eliminate the NAC terms from BO close-coupling equations was first discussed by Hobey and McLachlan26 for a single degree of freedom only and attempted by Smith for a diatomic molecule.27 M. Baer derived a general formulation28 on the diabatization of two adiabatic PESs for atom−molecule collisions along with ADT condition, where the ADT matrix elements were shown to be obtained by solving the differential equations for ADT angles along a two-dimensional contour, and the existence and uniqueness of the solution of those equations was predicted to be guaranteed by the validity of the curl condition.29,30 Received: November 25, 2012 Revised: March 9, 2013 Published: March 22, 2013 3475

dx.doi.org/10.1021/jp311597c | J. Phys. Chem. A 2013, 117, 3475−3495

The Journal of Physical Chemistry A

Article

Mead and Truhlar31 and Kendrick32 also discussed the removable and nonremovable components of NACTs, where the ADT can at best remove the longitudinal component of the derivative coupling.28,29 At the close proximity of a degeneracy, only the removable coupling is singular and the nonremovable coupling is insignificant.33 The ADT angle(s) can be calculated by integrating the derivative coupling at and around the CI point(s),29 where the closed line integrals of the derivative coupling will be multiples of π.34,35 On the other hand, when the nonremovable couplings are not negligible, the residue of such closed line integrals does not appear as multiples of π.36 Though the removable and nonremovable couplings can be separated by solving Poisson’s equation for the ADT angle (θ),37,38 the solution of Poisson’s equation is computationally expensive due to many possible definitions for the boundary conditions on such angles.39 Indeed, one can reduce the magnitude of nonremovable coupling by including a greater number of electronic states, which again increases the computational cost enormously. In any case, the requirements for diabatization through ADT around CI point(s) are the following: (a) the singularity in the derivative coupling must be transformed away and (b) the residual couplings must be negligible. Since the NACTs are usually very sharp functions of nuclear coordinates with singularities,40−42 one needs to perform a unitary transformation to obtain the diabatic representation of SE, where potential couplings among the electronic states are slowly varying functions of nuclear coordinates, and thereby, the dynamical calculations could be numerically accurate and stable. Such a transformation for a given sub-Hilbert space (i.e., an isolated group of states form a complete space) is guaranteed only when the vector fields created by the NACTs satisfy the so-called curl conditions29 at every point in the region of interest in order to eliminate the NACTs. It may be noted that this condition is not fulfilled at degenerate points (or CI points), and nevertheless, the NACTs can be eliminated if the line integrals are quantized forming a subspace. Sarkar et al.14−19 generalized the BO treatment for any three/four coupled electronic states in terms of electronic basis functions/ADT angles, where the explicit form of NACTs, their curl−divergence equations, curl condition, and diabatic PESs in terms of ADT angles were formulated. This approach provides a practical way to handle the NACTs and helps to construct continuous and single valued diabatic PESs. The spectroscopy of alkali trimers are important for exhibiting JT interaction, where their convenient oneelectron-valence configuration and small number of vibrational degrees of freedom (DOF) make them easy to analyze both by experimental measurements and theoretical calculations. In Li3 cluster, Sadygov and Yarkony predicted43 unusual CIs of C2v symmetry along with that of D3h symmetry. Koizumi and Bersuker modified44 the existing Pseudo Jahn−Teller (PJT) model Hamiltonian45 and predicted four CIs between the excited 22E′ states of Na3 cluster. We performed ab initio calculation for PESs and NACTs of 22E′ and 12A1′ states of Na3 cluster at the SA-CASSCF level and reported four CIs46−48 between the states of 22E′ symmetry. With these ab initio quantities, we had constructed diabatic surfaces and carried out dynamics to obtain the absorption spectrum of Na3 cluster. Since those diabatic PESs were not appearing as symmetric and smooth functions with respect to the normal modes (Qx and Qy) and the calculated absorption spectrum47 was away from the latest experimental one,49−54 we modify beyond BO theory

by employing MS and perform more accurate ab initio calculations with MRCI approach to construct diabatic surfaces as described below. Any ab initio quantum chemistry package assigns IREPs to the adiabatic PESs, but does not for the IREPs of the calculated NACTs. While constructing the diabatic PESs by incorporating multiple adiabatic PESs and NACTs, it is theoretically necessary that both quantities of adiabatic PESs and NACTs belong to appropriate IREPs. The only way to assign IREPs to the NACTs is to adapt them with MS as proposed by Longuet−Higgins,55 where, for degenerate states, we need to construct the generators of the character table of D3h(M) symmetry. We impose MS56 on NACTs by considering appropriate nodal planes for different IREPs and location of CIs. At present, MRCI approach has been employed to obtain both adiabatic PESs and NACTs to bring more accuracy. Ab initio calculations for adiabatic PESs and NACTs are performed as a function of the symmetric coordinate [ρ = ((Qx2 + Qy2))1/2 and ϕ = tan−1 (Qy/Qx)] to incorporate the inherent symmetry of the molecule. Both these theoretical developments and technical improvements are expected to help to construct single valued, continuous, smooth, and symmetric diabatic PESs, and the dynamical calculations on those surfaces could reproduce the pattern of the experimental two-photon ionization (TPI) spectrum of system B of Na3 cluster.

II. THEORY A. Beyond Born−Oppenheimer Theory for ThreeState Sub-Hilbert Space. We assume a three-state electronic sub-Hilbert space with CIs among the states anywhere in the nuclear CS. The BO expansion of the molecular wave function for this subspace of the Hilbert space and the total electron− nuclei Hamiltonian in the adiabatic representation are 3

Ψ(se , sn) =

∑ ψi(sn)ξi(se; sn) i=1

Ĥ = Tŝ n + Ĥ e(se , sn) Tŝ n = −

ℏ2 2

⎛∇

∑ ⎜⎜ i

2⎞

⎟ ⎟ ⎝ mi ⎠ sn , i

Ĥ e(se , sn)ξi(se ; sn) = ui(sn)ξi(se ; sn)

(1)

where T̂ sn is the nuclear kinetic energy (KE) operator and the expansion coefficient, ψi(sn), are actually the nuclear wave functions. The eigenfunctions [ξi(se;sn)] of the electronic Hamiltonian, Ĥ e(se, sn), are defined by the sets of nuclear (sn) and electronic (se) coordinates with nuclear coordinate dependent eigenvalues, ui(sn). From the time independent SE [Ĥ Ψ(se, sn) = EΨ(se, sn)] for the total electron−nuclear Hamiltonian and the BO expansion of the molecular wave function [eq 1], the elements of nonadiabatic coupling matrices of the first [τ⃗(1)] and second [τ(2)] kind are τij⃗(1) = ⟨ξi(se ; sn)|∇⃗ξj(se ; sn)⟩

(2a)

τij(2) = ⟨ξi(se ; sn)|∇2 ξj(se ; sn)⟩

(2b)

respectively. Since, for a given sub-Hilbert space, the matrices, τ⃗(1) and τ(2) are related as: 3476


The Journal of Physical Chemistry A τ (2) = τ ⃗(1)·τ ⃗(1) + ∇⃗τ ⃗(1)

Article

matrix (A) is constituted with orthogonal basis and satisfies the orthonormality condition [ATA = I], it belongs to a rotational group denoted as O(3), and such matrix with detA = +1 is called special orthogonal group, SO(3). Thus, any operation under the rotation matrix (A) should preserve the magnitudes of the functions and the angles between them, namely, adiabatic PESs. (b) Since the ADT angles (i) are not the DOF, at present case, they rather appear as the functions of the coordinates, ρ and ϕ (DOFs), and (ii) could rarely attend the exact magnitude of 0 or nπ as the numerical solution of ADT equation (eq 9), the Gimbal Lock in ADT matrix will not generally arise. Even if it arises, it will not hamper the construction of diabatic PESs with all required properties, namely, single-valued, continuous, symmetric, and smoothness of the surface because the numerical values of ADT angles are not arbitrary, but defined from the solution of ADT equation (eq 9). (c) The SO(3) group is nonabelian because the rotation of an angle around an axis followed by the rotation around another axis is different from the overall rotation in reverse order. In other words, the change of an ADT angle due to the integration along a coordinate (ρ or ϕ) followed by its change with the integration along another coordinate (ϕ or ρ) is different from the overall change of the ADT angle while performing the integration in reverse order. Such different paths of integration to obtain ADT angles lead to different ADT matrices, and each of the ADT matrices will, therefore, be related with the others through orthogonal transformations. The Appendix B and Supporting Information show the property of orthogonality between the two matrices due to the choice of two different paths. (d) The Lie algebra of SO(3) consists of a skew-symmetric matrix. The matrix equation for ADT (eq 7) involved with skew-symmetric matrix τ (NAC matrix) will generate a path-dependent solution (ADT angles) due to the choice of contour on the coordinates, ρ and ϕ. The corresponding ADT matrix will be related through orthogonal transformation with any other ADT matrix constructed by a different choice of contour. (e) The components of the ab initio calculated NACTs are substituted in ADT equation (eq 7), and the coupled stiff differential equations (eq 9) are solved by using backward differentiation formula (BDF) to calculate the ADT angles. Those angles are substituted in ADT matrix (eq A.1) to obtain the diabatic SE (eq 6), where the explicit form of the diabatic PES matrix elements in terms of ADT angles are shown in Appendix A. (f) When we transform the kinetically coupled adiabatic SE (eq 4) to the diabatic one (eq 6), it is necessary to make sure both the kinetic and potential energy terms in the adiabatic representation are consistent with the appropriate symmetry, i.e., totally symmetric under the operation of each symmetry element of MS group, D3h(M); otherwise, the diabatic SE will lead to erroneous solution. In other words, if we do not care about the IREPs of the NACTs calculated from ab initio packages, plug those terms into the ADT equation (eq 9) to obtain the ADT angles and construct the diabatic surfaces (see eq A.2); the corresponding diabatic surface elements47 could not show continuity, symmetry, and smoothness of their functional forms, and thereby, the theoretically calculated spectra was not in good agreement with the experimental one. However, because the ab initio quantum chemistry packages do not provide appropriate IREP adapted NACTs, we employ the socalled MS of the system to implement such description. The adiabatic PESs are the solution of electronic Hamiltonian by freezing the nuclear geometry at each grid point within the BO treatment. On the contrary, the NACTs are calculated by

(3)

we can arrive at the following compact form of kinetically coupled nuclear equations: −

ℏ2 ⃗ (∇ + τ ⃗)2 ψ + (U − E)ψ = 0 2m

(4)

with the following form of nonadiabatic coupling [τ⃗(≡ τ⃗(1))] and adiabatic PESs matrices: ⎛ 0 τ12⃗ τ13⃗ ⎞ ⎜ ⎟ 0 τ23 τ ⃗ = ⎜−τ12⃗ ⃗ ⎟ and Uij = uiδij ⎜⎜ ⎟⎟ 0⎠ ⃗ ⎝ −τ13⃗ −τ23

(5)

respectively. When the three states constitute a sub-Hilbert space, one can transform (ψ = Aψd) the adiabatic nuclear SE (eq 4) to the diabatic one: −

ℏ2 2 d ∇ ψ + (W − E)ψ d = 0, 2m

W = A†UA

(6)

with the condition: ∇⃗A + τ ⃗ A = 0

(7)

This equation is known as the adiabatic−diabatic transformation (ADT) condition.28 It is possible to find meaningful solution from eq 7 only when the chosen form of A matrix satisfies (a) the orthonormality condition at any point in CS and (b) periodicity with respect to a parameter. Since the model form of (3 × 3) A matrix has to be orthogonal, the orthonormality conditions demand the fulfillment of six relationships. Thereby, three independent variables commonly called ADT angles, viz., Euler like angles of rotation [θ12(sn), θ23(sn), and θ13(sn)], are the natural requirement to construct the three-state A matrix by taking the product of three (3 × 3) rotation matrices, A12(θ12), A23(θ23), and A13(θ13) in different order, where one of them is: A(θ12 , θ23 , θ13) = A12(θ12) ·A 23(θ23) ·A13(θ13)

(8)

When the antisymmetric form of τ matrix (eq 5) and the model form of A matrix (see eq A.1 in Appendix A) are substituted in eq 7, Top and Baer57 and Alijah and Baer35 arrived at the following equations for ADT angles: ∇⃗θ12 = −τ 12 ⃗ + tan θ23(τ ⃗ 13cos θ12 − τ ⃗ 23 sin θ12)

(9a)

∇⃗θ23 = −(τ 13 ⃗ sin θ12 + τ ⃗ 23 cos θ12)

(9b)

∇⃗θ13 = −

1 (τ 13 ⃗ cos θ12 − τ ⃗ 23sin θ12) cos θ23

(9c)

However, if we substitute the different columns of A matrix (see eq A.1 in Appendix A) in eq 2a, the explicit forms14 of τ matrix elements in terms of ADT angles are given by: τ12⃗ = −∇⃗θ12 − sin θ23∇⃗θ13

(10a)

τ23 ⃗ = sin θ12cos θ23∇⃗θ13 − cos θ12∇⃗θ23

(10b)

τ13⃗ = −cos θ12cos θ23∇⃗θ13 − sin θ12∇⃗θ23

(10c)

The beyond BO equations (eqs 4−10) have the following theoretical aspects and technical consequences: (a) As the ADT 3477



Article

Figure 1. Permutation of nuclear coordinates is denoted by the operation 123, and it is equivalent with the successive operation of C13, R2π/3, and p123. The operation C13 acts only on electronic and nuclear space coordinates rotating the vibrational displacements and electronic coordinates with a molecular fixed axis, R2π/3 is a bodily rotation of the molecule about the molecular fixed axis, and p123 successively permutes the spins of the nuclei.

Figure 2. (a) Symmetry elements of Na3 system described by point group D3h. (b) The effects of the symmetry operations in point group D3h.

Ô = Oâ Ô bOĉ

considering the nonrigidity of the nuclei and thereby, MS adaptation is a necessity to assign the IREPs of those terms for beyond Born−Oppenheimer treatment. B. Molecular Symmetry Adaptation to NACTs for D3h(M) Group. The true symmetry group for the complete molecular Hamiltonian is called molecular symmetry group, which, in general, contains R symmetry operations, Ô i, i = 1,2,...,R, including the identity operation E, all feasible permutations P of the spatial as well as spin coordinates of the equivalent nuclei, the inversion E* of all nuclear coordinates sk and electronic coordinates se, and the inversion-permutation P·E*. Each symmetry operation Ô of the MS group of a molecule transforms the vibronic modes, the three angles (Euler) to describe the orientations of a rigid body in a threedimensional Euclidean space and the nuclear spin, which enables us to express Ô as a product of three different operations:58

(11)

where Ô a produces change in vibronic DOF, Ô b changes the Euler angles, and Ô c permutes the nuclear spin. In contrast to the MS group, the symmetry operations of the point group of a molecule transform only the vibronic DOF. At this point, we take the example of sodium trimer that belongs to the point group D3h. If one locates the system-fixed (x,y,z) axes on a Na3 cluster such that the xy plane is the molecular plane (plane of symmetry, σh) and the z axis is the 3-fold symmetry axis (C3), the relationships between the MS group [D3h(M)] and the corresponding point group [D3h] operations are as follows58 (one representative element from each class of the MS group is shown):

3478

E = E ·R0·p0

(12a)

(123) = C31·R2π /3·p123

(12b)



Article

Table 1. Extended Character Table of Molecular Symmetry Group D3h(M) for a NAC Element of Different Nuclear Coordinates D3h(M)

E

(123) (132)

(12) (13) (23)

E*

(123)* (132)*

(12)* (13)* (23)*

A1′ A2′ E′ A″1 A″2 E″

1 1 2 1 1 2

1 1 −1 1 1 −1

1 −1 0 1 −1 0

1 1 2 −1 −1 −2

1 1 −1 −1 −1 1

1 −1 0 −1 1 0

coord

first comb

deriv

ϕ x,y

∂/∂ϕ ∂/∂x, ∂/∂y

ρ

∂/∂ρ

second comb

τρ

τϕ

τϕ

τx, τy τρ

τx,τy

Table 2. Extended Character Table of Molecular Symmetry Group D3h(M) with Different NAC Elements D3h(M)

E

(123) (132)

(12) (13) (23)

E*

(123)* (132)*

(12)* (13)* (23)*

A1′ A′2 E′ A″1 A″2 E″

1 1 2 1 1 2

1 1 −1 1 1 −1

1 −1 0 1 −1 0

1 1 2 −1 −1 −2

1 1 −1 −1 −1 1

1 −1 0 −1 1 0

coord

deriv

ϕ x,y

∂/∂ϕ ∂/∂x, ∂/∂y

ρ

∂/∂ρ

13 τ23 k τk

τϕ τx,τy τρ

τ12 k τϕ τx, τy τρ

Figure 3. Nodal patterns of the IREPs of molecular point group D3h.

(23) = C2x·Rπz·p23

(12c)

E* = σxy·Rπz·p0

(12d)

(132)* = S3−1·R zπ /3·p132

(12e)

(23)* = σxz·R yπ ·p23

(12f)

and C(3) 2 ) are in xy plane, each passing through one vertex of the triangle and bisecting the opposite side of the same. There (2) (3) are three vertical symmetry planes also, viz. σ(1) v , σv , and σv , each containing one C2 axis and is perpendicular to the xy plane or better to say the horizontal symmetry plane, σh (not shown in Figure 2a). As all the symmetry elements of the D3h point group can not be displayed in such a simple manner that is adapted in Figure 2a, we use Figure 2b to depict a modified version of the traditional projection59 of the same point group, where the effects of the symmetry operations are indicated without showing the symmetry elements. If one considers a point E above the xy plane representing the identity operation, the different symmetry operations acted separately upon that point would displace the same to different positions, and Figure 2b shows all such positions together indicating the appropriate symmetry operations in each case. There are two important theorems56 on MS, which relate the IREPs of different NACTs.

ϕ

where R is the member of the corresponding molecular rotation group and p denotes nuclear permutation. Figure 1 illustrates eq 12b schematically. With the three coordinate axes (x,y,z) being located on the Na3 system as described in Figure 1, we represent the symmetry elements of the system in Figure 2a, where the planar system is considered to be on the xy plane. The principal axis of symmetry, C3, is perpendicular to that plane (i.e., along the zaxis) passing through the centroid of the equilateral triangular (2) system and the other three rotational symmetry axes (C(1) 2 , C2 , 3479



Article

Figure 4. (a,b) Radial and torsional parts of a closed loop Lg around a CI located in (ρ, ϕ) plane for the first combination of IREPs of the NACTs. (c,d) Same for the second combination.

(i). Theorem 1. The first theorem states that when the IREP of a particular NACT, say τijl , is known for a specific symmetry adapted nuclear coordinate sl, the IREP of another one, τijk, between the same electronic states i,j for a different symmetry adapted nuclear coordinate, sk, can be found out by the following relationship: ⎛ ∂ ⎞ ⎛∂ ⎞ Γ(τkij) = Γ⎜ ⎟ × Γ⎜ ⎟ × Γ(τlij) ⎝ ∂s l ⎠ ⎝ ∂s k ⎠

respectively. As we perform the dynamics on the diabatic surfaces for a fixed positive Qs(=3.7 Å), it is necessary to consider only the smaller circles of Figure 3, while determining the nodal patterns of the NACTs. Since we need to select the correct combination from the two possible ones (see Table 1), we utilize the quantization rule60 (eq 15) of the contour integrals over the NACTs evaluated along a closed loop Lg of nuclear coordinate sn around a single conical intersection (CI):

(13)

∮ dsn·τ(sn|Lg) = ±π

where τkij = ⟨ξi(se ; sn)|

∂ (ξj(se ; sn))⟩ ∂sk

(15)

as shown in Figure 4. Such a rule was established and implemented60 for locating the number of CIs by analyzing the NAC elements over the contour integral, 0 to 2π. Thereby, the contour integral can be divided into three components,56 viz., one torsional (tor) and two radial (rad) lines that create the loop:

(14)

The possible combination of the IREPs of the NACTs can be determined by fixing each of them [τk, (k = ρ, ϕ, x, and y)] in succession to the totally symmetric representation (A′1). In other words, when we choose τρ = A′1, we get the other three as τϕ = A1″, τx = τy = E″ from Table 1 and eq 13. Though it seems that we should find four such combinations of the IREPs of the NACTs, only two possibilities survive (see the last two columns of Table 1) as theorem 1 is valid only for those Γ(∂/∂sl) that belong to 1-D IREP (i.e., for l = ρ, ϕ). On the contrary, when we try to find out the combination of the IREPs of the τks fixing τx or τy to A′1; the IREP for τy or τx turns into reducible representation contradicting the construction of the character table as shown in Tables 1 and 2. With the aid of Figure 2b and Table 1, we depict the nodal patterns of the NACTs belonging to all possible IREPs in Figure 3. The profiles of the NACTs are displayed as functions of the normal mode coordinates Qx (bending) and Qy (antisymmetric stretching), where the outer and inner concentric circles represent the patterns for negative and positive values of Qs (symmetric stretching), respectively. In other words, the smaller and the larger concentric circles in Figures 2b and 3 display points above and below the xy plane,

∮ dsn·τ(sn|Lg) = Irad,1 + Itor + Irad,2 = ±π

(16)

where Irad,1 =

∫0

Irad,2 =

ρ1

dρτρ(ρ , ϕ1|Lg ),

Itor =

∫ϕ

ϕ2

dϕτϕ(ρ1 , ϕ|Lg ),

1

∫ρ

0

dρτρ(ρ , ϕ2|Lg )

1

(17)

When we introduce the symmetry of the integrands (τρ/τϕ) as defined in Table 1, incorporate the sign of their magnitude according to the inner circle of Figure 3, and perform those integrations, it appears that (a) for the first combination, Irad,1 + Irad,2 = 0 (from Figure 4a) and Itor = 0 (from Figure 4b), which is incompatible with the quantization rule (eq 15); (b) the second combination gives nonzero residue for the contour integral (see Figure 4c,d) and thereby seems to be the only feasible one. 3480



Article

⎛ ⎞ ∂ Γ(τx23) = Γ⎜ ξ2(se ; sn)| (ξ3(se ; sn)) ⎟ ⎝ ⎠ ∂x

(ii). Theorem 2. To state the second theorem (see eq 20), we need to consider a loop-type sequence of N molecular states with the same spin multiplicity, say N doublet states, Da, Db, Dc,...,Dy,Dz=a, for the Na3 cluster. The corresponding NACTs b,c y,z=a are τa,b , and their product is: k , τk ,..., τk τka , b , c ,..., y , z = a = τka , bτkb , c ...τky , z = a

⎛∂ ⎞ = Γ(ξ2(se ; sn)) × Γ⎜ ⎟ × Γ(ξ3(se ; sn)) ⎝ ∂x ⎠

As Γ(∂/∂x) = E′ and Γ(τ23 x ) = E′ (see Table 2), it appears from eq 23 that

(18)

In terms of the IREPs, one can write:

Γ(ξ2(se ; sn)) × Γ(ξ3(se ; sn)) = A′1

(24a)

Similarly, for Γ(τ13 x ) = E′, we have

Γ(τka , b , c ,..., y , z = a) = Γ(τka , b) × Γ(τkb , c) × ... × Γ(τky , z = a)

Γ(ξ1(se ; sn)) × Γ(ξ3(se ; sn)) = A′1

= Γ(ξa(se ; sn))2 × Γ(ξb(se ; sn))2 × ... × Γ(ξy(se ; sn))2

(24b)

It is clear from eqs 24a and 24b that ξ1(se;sn), ξ2(se;sn), and ξ3(se;sn) must belong to 1-D IREP, which means

⎛ ∂ ⎞N × Γ⎜ ⎟ ⎝ ∂s k ⎠

(a). For 1-D IREPs.

(23)

(19)

Γ(ξ1(se ; sn))2 × Γ(ξ2(se ; sn))2 × Γ(ξ3(se ; sn))2 = A′1

56

⎧ ⎛ ⎞ ⎪ Γ⎜ ∂ ⎟ when N is odd ⎪ Γ(τka , b , c ,..., y , z = a) = ⎨ ⎝ ∂s k ⎠ ⎪ ⎪ A′ when N is even ⎩ 1

When we use eq 19 for the present case and substitute eq 25, we obtain Γ(τx1,2,3,1) = Γ(ξ1(se ; sn))2 × Γ(ξ2(se ; sn))2 × Γ(ξ3(se ; sn))2 ⎛ ∂ ⎞3 ⎛ ∂ ⎞3 ⎛ ∂ ⎞3 × Γ⎜ ⎟ = A′1 × Γ⎜ ⎟ = Γ⎜ ⎟ = (E′)3 ⎝ ∂x ⎠ ⎝ ∂x ⎠ ⎝ ∂x ⎠

(20)

Γ(τx1,2,3,1) = Γ(τx12) × Γ(τx23)2 = Γ(τx12) × (E′)2

(27)

Γ(τ12 x )

By equating eqs 26 and 27, one obtains = E′, and In Table 2, we present the IREPs for all ) = E′. similarly, Γ(τ12 y kinds of NACTs in a compact manner.

III. PESs AND NACTs, LOCATION OF CIs, AND MS ADAPTED NACTs A. Ab Initio Calculation for PESs and NACTs. The vibrational modes Qs, Qx, and Qy are referred to the symmetric stretching, bending, and antisymmetric stretching mode, respectively. The breathing mode Q s preserves the D3h symmetry, whereas Qx and Qy distort the D3h equilateral triangle into a C2v isosceles and a Cs scalene triangle, respectively. If the antisymmetric stretching coordinate, Qy, is set to zero, the states under D3h symmetry split into the corresponding C2v states, i.e., the potential energy curves of 22E′ and 12A′1 symmetry give rise to 22A1, 22B2 and 32A1 states at various nonzero values of Qx. We have employed a (3s,3p) Gaussian basis set, known as triple-ζ or TZV basis in our ab initio calculations, where such bases are being augmented with a single set of d functions (ζd = 0.1) to represent the electron correlation effects and two diffuse s functions (ζ4s = 0.008736; ζ5s = 0.00345)45 to accommodate the excited state description. The excited state PESs (22E′ and 12A′1) and NACTs calculated at SA-CASSCF level were already presented46−48 as functions of the normal modes Qx and Qy for three fixed Qs, namely, 3.4, 3.7, and 4.0 Å, respectively. In this article, we repeat those calculations with the more accurate methodology, namely, MRCI, by using the MOLPRO quantum chemistry package61 for a fixed Qs = 3.7 Å. Since 22E′ and 12A′1 states of Na3 cluster are actually the fourth, fifth, and sixth excited states of the system, it is necessary to show that they form a sub-Hilbert space within the range of CS of our investigation. MRCI ab initio calculations are performed with the first six states, where those reference states are calculated at first by the MCSCF approach involving the first 20 electronic

(21)

To be specific, let the excited electronic states 22E′ and 12A′1 of Na3 system be denoted as 1 ,2, and 3, respectively, where multiple CIs are predicted46−48 between states 1 and 2. Thus, 23 τ13 k and τk should interchange their respective signs while crossing each CI, and they should have the same IREP. According to eq 19, for a three-state sub-Hilbert space like that created by the states 22E′ and 12A1′ in Na3 system, one can write56 the following product for the IREPs of the NACTs considering a loop -type sequence of those three states: Γ(τk1,2,3,1) = Γ(τk12) × Γ(τk23) × Γ(τk31) = Γ(τk12) × Γ(τk23)2

(26)

as from Table 2, Γ(∂/∂x) = E′. From Table 2 and eq 22

While implementing the second theorem (eq 20) on the NACTs, we intend to discuss the following property56 of the wave functions of a system near a CI. If there is a CI between two electronic states i and (i + 1), the corresponding electronic eigenfunctions ξi(se;sn) and ξi+1(se;sn) may change their signs along the nuclear coordinate sn at the time of crossing the CI. When all the CIs in CS are connected with a symmetry related set of coordinates and there is an interchange of sign in the wave functions for one of the CIs, similar switching will occur for all the other CIs, too,56 and thereby, the wave functions ξi(se;sn) and ξi+1(se;sn) will have the same IREP. As a consequence, the NACTs τi,jk and τ(i+1),j , i < j, should also k interchange their signs (see eq 14) with the same IREP, i.e., Γ(τki , j) = Γ(τk(i + 1), j)

(25)

(22)

where Γ(τijk) = Γ(τjik) for any i,j, and the so-called crossing rule (eq 21) between the states i and (i + 1) are incorporated. 23 If we consider that Γ(τ23 ρ ) = A″ 1 and Γ(τϕ ) = A′1 (see Table 12 12 2), one obtains Γ(τρ ) = A2″ and Γ(τϕ ) = A2′ from eq 22, where Γ(τ1,2,3,1 ) = Γ(∂/∂ρ) = A2″ and Γ(τ1,2,3,1 ) = Γ(∂/∂ϕ) = A2′ ρ ϕ according to eq 20. (b). For 2-D IREPs. Since eq 20 is only valid to find out the 1D IREPs, we establish similar working equations for 2-D IREP cases as given below: 3481



Article

Figure 5. Contour plot of the PESs, where the surfaces at each point are scaled with the minimum energy of the lower 22E′ state (−482.5668 hartree). (a,b) Lower and upper 22E′ state; (c) 12A′1 state.

components of NACTs are projected in normal mode coordinates x(Qx) and y(Qy) as shown in Figure 6. At this point, we wish to show that the CI at ρ = 0 could be a consequence of the collapse of three CIs. Figure 7 displays the following: (a) at relatively higher ρ value (1.3 × 10−4 Å), the magnitude of τ12 ρ appears zero after the π/3 interval of the ϕ coordinate; (b) on the contrary, as the ρ coordinate tends to zero (0.1 × 10−4 Å), such magnitude reduces to a lower value but not zero on the same interval of ϕ. B. Longuet−Higgins’ Phase Change Theorem to Locate CI(s). In order to locate the CIs in the CS, we also calculate the Longuet−Higgins’ (LH)20,21 sign change by following the numerical magnitude of the component of the dominant configuration around a close loop of one or multiple CI(s).24 Figure 8 depicts (I) the four loops that successively encircle one or multiple CI(s) and (II) their corresponding sign changes. The figure shows that, while adiabatically transporting along a close path encircling odd or even number of CIs, the MRCI wave functions do or do not change its sign as predicted by Longuet−Higgins’ (LH) theorem. C. Sign Correction on NACTs Employing MS. Since it is well-known that the existing quantum chemistry packages still

states to reach the closest possible variational minimum with respect to the exact energies. For evaluating the excited state PESs of symmetry 22E′ and 12A1′ , 15 orbitals are considered to be closed shell and the three unpaired electrons are distributed over 12 active orbitals generating 572 numbers of configuration state functions (CSFs). It could be noted here that the energy of the lowest 22E′ state obtained from MRCI level calculation is lower by ∼1538.95 cm−1 than the corresponding energy calculated by CASSCF approach. Figure 5 presents the contour plots of PESs for 22E′ and 12A′1, where the embedded energies are the relative energy values given in hartree unit. The energies are scaled by subtracting the minimum energy (−482.5668 hartree) of the lower 22E′ state. The NACTs among 22E′ and 12A1′ states are calculated by using the so-called numerical method (DDR) in MOLPRO. All these MRCI calculations for the adiabatic PESs and NACTs were performed over the range ρ = 0.0 to 0.9 Å and ϕ = 0 to 2π with a grid size of 150 × 150. As reported earlier, present calculations also show that there are four CIs between the 22E′ states of Na3 system, one of them being at (Qx,Qy) = (0,0) [(ρ, ϕ) = (0,0)] and the other three at ρ ≈ 0.03 Å, ϕ = π/3, π, and 5π/3, where ρ = (Qx2 + Qy2)1/2 and ϕ = tan−1 (Qy/Qx). The functional forms of the ρ and ϕ 3482



Article

Figure 6. Functional forms of the components of nonadiabatic coupling (NAC) matrix elements as functions of Qx and Qy. (a,d) The outermost contour line (dotted) starts with 0 and the inner ones (solid) rise up with an increment of 6 Å−1. (b,c,e,f) The outermost contour line (dotted) starts with 0 and the inner ones (solid) rise up with an increment of 0.2 Å−1.

crossing, which is a basic axiom (see eq 21) for establishing the second theorem of MS. Figure 9 presents τijk (k = ρ, ϕ) vs ϕ for a constant ρ = 0.03 Å, where the upper panels (a,b) display the numerical magnitude of τijρ and τijϕ as obtained from the MOLPRO package, and the lower panels (c,d) depict the IREP adapted (dotted line) ones with appropriate crossing due to CIs

can not provide appropriate IREPs of the NACTs, we have prepared MS-adapted NACTs through the following two steps: (a) a broad overall sign correction is performed for each of the NACTs according to the nodal patterns of the respective IREPs 23 to which they belong (see Figure 3); (b) the signs of τ13 k and τk [k = ρ,ϕ] are interchanged at the point of CI to incorporate 3483



Article

achieve the required level of convergence. It is important to note that one can integrate two sets of coupled (ρ and ϕ) differential equations (eq 9) with infinitely different chosen paths, where each path will produce a set of ρ−ϕ dependent ADT angles and such angles obtained from any choice of path are expected to show gauge invariance. Figure 10 describes two such rectangular paths: S.I., The differential equations for the ϕ grid are integrated with positive increment from ϕ = 0 to 2π for each positive step integration of the differential equations for the ρ grid from 0 to ρmax; S.II., The differential equations for the ρ grid are integrated with positive increment from ρ = 0 to ρmax for each positive step integration of the differential equations for the ϕ grid from 0 to 2π. We display the ADT angles [θij(ρ,ϕ)] as functions of ρ and ϕ in panels a−c of Figures 11 and 12 obtained by using S.I and S.II schemes, respectively. The residues of the ADT angles [θ̅ikj = ∫ 02πθ ij(ρk,ϕ)dϕ] and their cumulative residues [θ̅ij,l = − θ̅kij|] are presented as a function of ρ in panels ∑kl = 1|θ̅k+1 ij d−f of Figures 11 and 12. The numerical magnitude of θ̅kij essentially counts the number of CIs between a pair of states and exhibits the gauge invariance for two paths of integration, S.I and S.II. When a system has multiple CIs, the contribution for each CI on residue (θ̅kij) could be positive or negative, and thereby, the cumulative residue (θ̅ij,l) could be a better way of demonstrating the total change of the residue. Panel d of Figures 11 and 12 reveals that the cumulative residue (θ̅ij,l) for the ADT angle θ12 become close to 4π at asymptotic ρ indicating the presence of four CIs between the 22E′ states, where the one to one correspondence between the panels d−f of Figures 11 and 12 for each cumulative residue (θ̅ij,l) depict the gauge invariance on the paths of integration, S.I−S.II. The residue in the panels d, e, or f of the Figures 11 and 12 is close to either 4π or π/2 implying the existence of gauge invariance. Similarly, any other diagonal path for contour integration is expected to show the gauge invariance of ADT angles with the same residue or cumulative residue. When we substitute those ADT angles calculated by using a specific path to construct ADT (eq A.1) and diabatic PES (eq A.2) matrices, such matrices would be related through an orthogonal transformation with the corresponding matrices obtained from any

Figure 7. NAC element τ12 ρ as function of ϕ for different values of ρ(×10−4): (a) 0.1, (b) 0.3, (c) 0.5, (d) 0.7, (e) 0.9, (f) 1.1, and (g) 1.3 Å.

(solid line). In summary, MS determines only the “+” or “−” sign of the NACTs where the magnitudes of the NACTs are obtained from ab initio calculation. Though the NACTs should strictly follow the inherent symmetry property of the system even with large magnitude at and around the CI point(s), in practice, the ab initio calculation provides only approximate values around those critical point(s) no matter how much sophisticated and accurate electronic structure calculations are performed. However, it is important to note that the singularity of the NACTs and their sign in the CS are much more important to describe the system accurately through diabatic surfaces rather than their absolute magnitude.

IV. ADT ANGLES (θIJ) AND THEIR RESIDUES (θ̅KIJ) We substitute the MS adapted NACTs in eq 9 and solve over a 2D grid of geometries in polar coordinates (ρ, ϕ) ranging from ρ = 0 to 0.9 Å and ϕ = 0 to 2π spanning 512 × 512 grid points to obtain the ADT angles. Since such equations are stiff differential ones, backward differentiation formula is employed with appropriate relative and absolute error tolerances to

Figure 8. (a) Close paths (circle) encircling one or multiple conical intersections. (b) ci, the coefficient of the dominant configuration in MRCI wave function of the ground state. 3484



Article

Figure 9. Upper panels (a,b) show τij vs ϕ for a fixed ρ = 0.03 Å before MS adaptation, and lower panels (c,d) represent the same after MS adaptation.

Figure 10. Different 2D rectangular contour paths along which the ADT angles are solved.

equation14 in terms of ADT angles for each NACT was obtained by using eqs 10 and 28, and the divergences17 of τ⃗ijs were derived by employing eq 10 as shown in Appendix C. As ∇⃗ θij and, in general, ▽2θij are nonzero at and around a CI in nuclear CS, it appears from eq C.2 that the divergence of the vector field τ⃗ij should also be nonzero for any arbitrary value of θij, which means that the vector field should have a nonzero curl too. If two electronic states i and (i + 1) have CI(s) between them, the NACTs of type τ⃗i(i+1)(s) could have singularity (pole) at the CI point(s) and they (the NACTs) decay as 1/r, r being the distance from the CI.38 Such a vector field can be resolved into two components, viz., longitudinal and transverse,39,63 where the curl of the former is, by definition, zero but that of the latter may or may not. Though the Abelian (longitudinal) and non-Abelian (transverse) magnitude of the curl equations are the key issue for the formulation of extended Born−Oppenheimer (EBO) equation,12,14 we need to explore the validity of the curl condition14−18 for a given sub-Hilbert space to construct the diabatic potential energy matrix. It may be noted that curl

other path, but the calculated observables for all cases are predicted to be same.

V. CURL−DIVERGENCE EQUATIONS A curl condition29 for each NACT, τ⃗ij, has been derived and proved to exist for an isolated group of states (i.e., a sub-Hilbert space) by considering the analyticity of the ADT matrix A for a pair of nuclear DOF ∂ ij ∂ ij τq − τp = (τqτp)ij − (τpτq)ij , ∂p ∂q τijq = ⟨ξi|∇q ξj⟩

τijp = ⟨ξi|∇p ξj⟩, (28)

where the curl due to vector product of NACTs and analyticity of ADT matrix are given by Cij = (τqτp)ij − (τpτq)ij and Zij = (∂/∂p)τijq − (∂/∂q)τijp, respectively, with ▽p = ∂/∂p and ▽q = ∂/∂q. The Cartesian coordinates p and q denote nuclear DOFs. For a two-state sub-Hilbert space, the curl−divergence equations are formulated62 for numerical calculations to explore the existence of sub-Hilbert space with a group of states. However, for the three-state case, the explicit form of curl 3485



Article

Figure 11. Left panels: The form of the ADT angles solved on a 2D contour integral by using S.I scheme. Right panels: The residues of the l k+1 k corresponding ADT angles [θ̅kij = ∫ 2π 0 θij(ρk,ϕ)dϕ] and the cumulative residues [∑k = 1|θ̅ij − θ̅ij|] as a function of ρ.

existence of the three-state sub-Hilbert space consisting of 22E′ and 12A1′ states of Na3 system.

condition is independent of Abelian or non-Abelian nature of the curl of the NACTs. We display the root-mean-square deviation (rmsd) between the mthematical (Zij) and ADT (Cij) curl employing the equation: χcurlτij(ρk) = (1/N)(∑Nl (Cij(ρk,ϕl) − Zij(ρk,ϕl))2)1/2 i,j = 1,2,3 to explore the numerical validity of the curl condition (eq 28). Since it is hard to visualize the differences between the two types of curls as functions of ρ and ϕ, we have scaled both the quantities, i.e., the mathematical (Zij) and the ADT (Cij) curls by the maximum magnitude of either one of them and depict their rmsd values. We calculate both the mathematical and ADT curls for each NACT, evaluate the deviation between the two kinds of curls and then add the square of such deviations along the ϕ coordinate to obtain the rmsds as functions of ρ. Since the calculated rmsd values for curl τ12 appear as absolutely zero, the other two rmsds for curl τ13 and curl τ23 are presented in Figure 13a,b for schemes S.I and S.II, respectively, where their negligibly small magnitude describe the validity of the curl condition and thereby show the

VI. ADT AND DIABATIC POTENTIAL MATRIX While solving the JT CI model for two electronic states by HLH, the analytic expression for the diagonal elements of the ADT matrix shows sign change along a contour around the CI point/seam, and such scheme was utilized12,13 and implemented60 for multistate model systems also. When we incorporate those ADT angles in eq A.1, diagonal element(s) of the ADT matrix, A11/A22 shows six sign changes in Figure 14a,c, and A33 displays no sign change in Figure 14b,d as function of ϕ at asymptotic ρ. Though six sign changes in the A11/A22 element of the A matrix could be due to six CIs between the E′ states, the residue/cumulative residue (see panels d−f in Figures 11 and 12) of the ADT angle θ12 indicates four CIs in the CS (see Figures 11d and 12d). One possibility64 could be due to the collapse of three central CIs on a single point at ρ = 0 (see Figure 7). Since the inner CIs arise 3486



Article

Figure 12. Left panels: The form of the ADT angles solved on a 2D contour integral by using S.II scheme. Right panels: The residues of the l k+1 k corresponding ADT angles [θ̅kij = ∫ 2π 0 θij(ρk,ϕ)dϕ] and the cumulative residues [∑k = 1|θ̅ij − θ̅ij|] as a function of ρ.

due to the symmetry of the molecule, D3h, there is no option left except that the three CIs at small ρ and different ϕ collapse at a single point. It may be noted even that Longuet−Higgins’ sign change technique can not be implemented here to locate the number of CIs at the central point (ρ = 0; ϕ = 0). Figures 15 and 16 depict the functional nature of the diabatic matrix elements calculated by substituting the θij in eq A.2, where those ADT angles are obtained by integrating eq 9 for two paths of integration, S.I and S.II, respectively. The diabatic surfaces appear not only continuous and single valued but also smooth and symmetric, particularly, by considering their crossing with each other along a line in the CS (see panels b−c and d−e in Figure 15). There are altogether six different ways to take the product of the three rotation matrices [A12(θ12), A13(θ13), and A23(θ23)] to obtain the ADT matrix (A):17,35 A = Pn{A12(θ12) ·A 23(θ23) ·A13(θ13)}

n = 1, ...N!

where Pn is the nth permutation between two rotation matrix and N is the size of the Hilbert space. Each ADT matrix provides different sets of differential equations for ADT angles (eq 9), NACTs (eq 10), and their curl equations [eq C.1]. However, if we work with any one of the above six sets of A matrix, the calculated numerical magnitudes of A matrix elements remain same. To illustrate this, we take another product order of the rotation matrices: A(θ12 , θ13 , θ23) = A12(θ12) ·A13(θ13) ·A 23(θ23)

(30)

With this model form of A matrix (see eq A.3 in Appendix A) and the antisymmetric form of τ matrix (eq 5), we get a different set of differential equations for the ADT angles employing the ADT condition (eq 7) as before: ∇⃗θ12 = −τ 12 ⃗ − tan θ13(τ 13 ⃗ sin θ12 + τ ⃗ 23 cos θ12)

(29) 3487

(31a)



Article

leading to the following analytic expressions of NACTs in terms of the ADT angles:

1 (τ 13 ⃗ sin θ12 + τ ⃗ 23 cos θ12) cos θ13

∇⃗θ13 = −τ 13 ⃗ cos θ12 + τ ⃗ 23 sin θ12

(32a)

τ23 ⃗ = −cos θ12 cos θ13∇⃗θ23 + sin θ12∇⃗θ13

(32b)

τ13⃗ = −cos θ13 sin θ12∇⃗θ23 − cos θ12∇⃗θ13

(32c)

Two sets of ADT angles for the same nuclear configuration are obtained by plugging the ab initio calculated MS adapted NACTs τ⃗12, τ⃗23, and τ⃗13 in succession into eqs 9 and 31. When we calculate the A matrix elements with these two different sets of ADT angles, the numerical magnitudes for the two sets, say [Aij(ρk, ϕl)]1 and [Aij(ρk, ϕl)]2, match exactly. The rmsds between the elements of A matrix for the two cases are presented in Table 3, where the rmsd is defined as χAij(ρk) = (1/N)(∑Nl ([Aij(ρk, ϕl)]1 − [Aij(ρk, ϕl)]2)2)1/2, with i,j = 1,2,3 for eight different ρ values, ρ = 0, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, and 0.14 Å. Such rmsds clearly indicate the closeness between the two sets of ADT matrix elements, [Aij]1 and [Aij]2. Therefore, it really does not matter which product order of the rotation matrices has been considered to construct the A matrix and the same is true for the diabatic surfaces [Wij]. Since by definition the W matrix is generated as W = A†UA, U being the adiabatic PES, the numerical magnitudes of the W matrix elements (see Table 4) for the above two sets of numerically different ADT angles exactly match with each other even though the analytic forms of the W matrix elements in terms of the ADT angles appear different (see Appendix A).

Figure 13. (a) The rmsds between mathematical (Zij) and ADT (Cij) curls for τ13 and τ23 by using S.I scheme and (b) the same obtained from S.II scheme.

∇⃗θ23 = −

τ12⃗ = −∇⃗θ12 + sin θ13∇⃗θ23

(31b) (31c)

Figure 14. The (1,1)/(2,2) and (3,3) elements of adiabatic-to-diabatic transformation (ADT) matrix are shown. The upper panels (a, b) show the ADT matrix elements obtained from the ADT angles solved on a 2D contour integral by using S.I scheme and the lower panels (c,d) display the same obtained by employing S.II scheme. 3488



Article

Figure 15. Contour plots of diabatic PESs (eV) and the coupling elements (eV) calculated by incorporating ADT angles into eq A.2, where those angles are obtained by solving eq 9 on a 2D contour integral by using S.I scheme.

VII. NUCLEAR DYNAMICS AND ABSORPTION SPECTRA Since it is well-known that the theoretically calculated spectrum65,66 for system B of Na3 cluster agrees quite well with experimental one for Qs = 3.7 Å, nuclear dynamics has been performed on the diabatic surfaces over a 2D grid with the normal mode coordinates, Qx and Qy. The time dependent nuclear wave function obtained by FFT−Lanczos method is used to calculate the autocorrelation function (C(t)). Fourier transform of C(t) gives absorption spectra

where

∫−∞ C(t ) exp(iωt )dt

(34)

⎛t ⎞ ⎛t ⎞ = ψ *⎜ ⎟ ψ ⎜ ⎟ ⎝2⎠ ⎝2⎠

(35)

Equation 35 is more accurate, computationally faster, and convenient to implement than the previous one (eq 34), but it is valid only if the initial wave function is real and the Hamiltonian is symmetric. The general form of the Hamiltonian, Ĥ (Qx, Qy) = T̂ n·1 + W, is defined with the kinetic energy operator, T̂ n = −(ℏ/2)((∂/ ∂Qx2) + (∂/∂Qy2)) and the diabatic potential matrix, W(Qx, Qy)

∞

I(ω) ∝ ω

C(t ) = ⟨ψ (t )|ψ (0)⟩

(33) 3489



Article

Figure 16. Contour plots of diabatic PESs (eV) and the coupling elements (eV) calculated by incorporating ADT angles into eq A.2, where those angles are obtained by solving eq 9 on a 2D contour integral by using S.II scheme.

Table 3. rmsd Values of the Two Different Sets of A Matrix Elements ρ (Å) 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14

χA11(ρk) 0.2 1.3 1.1 1.2 1.2 1.3 1.1 1.4

× × × × × × × ×

10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8

χA22(ρk) 0.0 1.1 1.1 1.3 1.2 1.3 1.2 1.5

× × × × × × ×

10−8 10−8 10−8 10−8 10−8 10−8 10−8

χA33(ρk) 0.0 0.0 0.0 0.3 0.7 0.7 0.9 1.0

× × × × ×

10−8 10−8 10−8 10−8 10−8

3490

χA12(ρk) 0.2 0.7 1.2 1.2 1.3 1.5 1.6 1.6

× × × × × × × ×

10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8

χA13(ρk) 0.0 0.4 0.6 0.6 0.8 0.8 1.0 1.2

× × × × × × ×

10−8 10−8 10−8 10−8 10−8 10−8 10−8

χA23(ρk) 0.0 0.3 0.5 0.8 0.9 0.9 1.1 1.4

× × × × × × ×

10−8 10−8 10−8 10−8 10−8 10−8 10−8



Article

Table 4. rmsd Values of the Two Different Sets of W Matrix Elements ρ (Å) 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14

χW11(ρk) 0.0 0.0 0.0 0.3 0.2 0.0 0.4 0.2

−8

× 10 × 10−8 −8

× 10 × 10−8

χW22(ρk) 0.0 0.0 0.0 0.2 0.0 0.2 0.2 0.2

−8

× 10

× 10−8 × 10−8 × 10−8

χW33(ρk)

χW12(ρk)

χW13(ρk)

0.0 0.0 0.0 0.0 0.0 0.3 × 10−8 0.3 × 10−8 0.5 × 10−8

0.0 0.0 0.0 0.0 0.3 × 10−8 0.3 × 10−8 0.0 0.2 × 10−8

0.0 0.0 0.2 × 10−8 0.0 0.0 0.0 0.2 × 10−8 0.2 × 10−8

χW23(ρk) 0.0 0.0 0.3 0.3 0.2 0.2 0.2 0.3

× × × × × ×

10−8 10−8 10−8 10−8 10−8 10−8

(see Figures 15 and 16). We consider a Gaussian wavepacket (GWP) as the time-independent solution of the Hamiltonian [Ĥ (Qx, Qy) = T̂ n + (1/2)(Qx2 + Qy2)] and perform the adiabatic−diabatic transformation with the ADT matrix to obtain the diabatic wave functions with different initial positions of the GWP in the column matrix. As for example, the diabatic wave function starting with the initial wave function located on the ground state is given by ψdia = Aψadia with ψadia = (GWP,0,0). We carry out the dynamics thrice locating the initial wave function on each of the three states, 22E′ (states 1 and 2) and 12A′1 (state 3), and then convolute the three results considering only a small contribution from the state 12A1′ to achieve the final spectrum. The crossings among the diabatic population profiles as a function of time (see Supporting Information) indicate the presence of CIs between the states. The diabatic densities spread symmetrically over the space with time and clearly reflect (see Supporting Information) the symmetry of the diabatic surfaces. While performing the dynamics for each of the three cases, we have expanded the wave functions by 128 × 128 plane wave basis functions and such number of basis functions are enough to reach the convergence with respect to the 512 × 512 one. Figure 17 presents the direct comparison of that calculated absorption spectrum with the experimental outcome obtained by Ernst and Golonzka,54 where they measured the spectrum through cw laser excitation followed by laser ionization and mass selective detection in a molecular beam experiment for a limited wavelength region. The absorption spectra in Figure 17b,c are obtained by performing the dynamics on those diabatic surfaces constructed by the ADT angles obtained by integrating eq 9 for two paths of integration, S.I and S.II, respectively. The similarities of the two spectra in the peak positions as well as in the intensity pattern are quite encouraging.

VIII. SUMMARY If few electronic states of a molecular system constitute a subHilbert space over a range of CS, theoretical validity of ADT is guaranteed by curl conditions,35 and thereby, one can construct the diabatic Hamiltonian within the same range. We have calculated adiabatic PESs and NACTs for the states 22E′ and 12A1′ of Na3 system at MRCI level using quantum chemistry package MOLPRO, adapted MS on those NACTs, and solved the differential equations for the ADT angles (θij)57 by invoking a stiff equation solver technique [backward differentiation formula (BDF)] over a 2D grid on normal modes, Qx and Qy (≡ ρ, ϕ) for a fixed Qs. With those ADT angles and the adiabatic PESs (obtained from ab initio calculation), we have constructed the diabatic potential energy matrix by employing the analytic expressions of its elements in terms of ADT angles [see Appendix A].

Figure 17. Direct comparison of the absorption spectra of the present theoretical calculation with the experimental TPI spectrum, measured by Ernst and Golonzka54 (a). Theoretical absorption spectra are obtained by performing dynamics on those diabatic PESs calculated by S.I (b) and S.II (c) schemes, respectively.

We numerically demonstrate the collapse of three conical intersections at the central point of those normal modes (Qx = 0; Qy = 0), where the diagonal elements of ADT matrix (A11/ 3491



Article

A22) show the expected number of sign change due to the CI points. Previously, this feature was absent47 with less accurate NACTs data calculated by SA-CASSCF methodology. We have also shown numerically the gauge invariance for the various path dependent solution (ADT angles) of ADT equation, which leads to different ADT matrices [special orthogonal group, SO(3)], and those matrices are related with each other through orthogonal transformations. This became possible due to the incorporation of appropriate MS to the ab initio calculated NACTs (known as sign correction). There are three aspects that make the presently calculated diabatic PES single valued, continuous, smooth, and symmetric: (i) incorporation of inherent symmetry by performing ab initio calculation in symmetric coordinate [ρ = (Qx2 + Qy2)1/2 and ϕ = tan−1 (Qy/ Qx)] system; (ii) accurate ab initio calculated data (particularly the NACTs) by using MRCI methodology; (iii) appropriate IREP adaptation to the ab initio calculated NACTs. Nuclear dynamics has been carried out on those diabatic surfaces with an aim to reproduce the experimentally measured TPI spectrum of system B of Na3 cluster. In this context, it is important to note that Paul et al.47 calculated the theoretical spectrum by performing the dynamics on the diabatic surfaces constructed with ab initio (SA-CASSCF methodology) adiabatic PESs and without MS-adapted NACTs. Such calculations were carried out for two different situations and observed results are the spectrum for the pseudo Jahn−Teller (PJT) condition grossly matched with the experimental outcomes,49,54 but that for the full Hamiltonian showed very little agreement even with the old experimental data.49 It has been expected that the reason might be the less accurate ab initio (SA-CASSCF) calculated NACT and their improper IREP. The present work, which performs dynamics on the full diabatic Hamiltonian constructed with ab initio (MRCI approach) adiabatic PESs and molecular symmetry adapted NACTs, reproduces the pattern of latest experimental TPI spectrum by Ernst and Golonzka54 reasonably well. Therefore, the necessity of MS-adapted NACTs to construct diabatic surfaces is quite justified.

where the diabatic potential matrix elements in terms of ADT angles are as follows:

APPENDIX A: DIABATIC POTENTIAL ENERGY MATRIX The functional form of the A matrix by choosing the following order of the product of the three rotation matrices as shown in eq 8 is given by:

(A.2f)

W11 = u 2 cos 2 θ23 sin 2 θ12 + u3(cos θ12 sin θ13 + cos θ13 sin θ12 sin θ23)2 + u1(cos θ12 cos θ13 − sin θ12 sin θ13 sin θ23)2 (A.2a)

W22 = u 2 cos2 θ12 cos2 θ23 + u3( −sin θ12 sin θ13 + cos θ12 cos θ13 sin θ23)2 + u1( −cos θ13 sin θ12 − cos θ12 sin θ13 sin θ23)2 (A.2b)

W33 = u3 cos2 θ13 cos2 θ23 + u1 cos 2 θ23 sin 2 θ13 + u 2 sin 2 θ23 W12 = W21 = u 2 cos θ12 cos2 θ23 sin θ12 + u3[(− sin θ12 sin θ13 + cos θ12 cos θ13 sin θ23) × (cos θ12 sin θ13 + cos θ13 sin θ12 sin θ23)] − u1[(cos θ13 sin θ12 + cos θ12 sin θ13 sin θ23) × (cos θ12 cos θ13 − sin θ12 sin θ13 sin θ23)]

(A.2d)

W13 = W31 = cos θ23[−(u1 − u3)cos θ12 cos θ13 sin θ13 + sin θ12(−u 2 + u3 cos2 θ13 + u1 sin 2 θ13)sin θ23] (A.2e)

W23 = W32 = −cos θ23[−(u1 − u3)cos θ13 sin θ12 sin θ13

■

− cos θ12(u 2 − u1 sin 2 θ13 − u3 cos2 θ13)sin θ23]

On the other hand, when the A matrix is constructed with the following order as introduced in eq 30: A(θ12 , θ13 , θ23) = A12(θ12)·A13(θ13)·A 23(θ23) ⎛ cos θ12 sin θ12 0 ⎞⎛ cos θ13 0 sin θ13 ⎞ ⎟ ⎜ ⎟⎜ = ⎜− sin θ12 cos θ12 0 ⎟⎜ 0 1 0 ⎟ ⎜ ⎟⎜ ⎟ ⎝ 0 0 1 ⎠⎝− sin θ13 0 cos θ13 ⎠ ⎛1 ⎞ 0 0 ⎜ ⎟ ⎜ 0 cos θ23 sin θ23 ⎟ ⎜ ⎟ ⎝ 0 − sin θ23 cos θ23 ⎠

A(θ12 , θ23 , θ13) = A12(θ12)·A 23(θ23)·A13(θ13) ⎛ cos θ12 ⎜ = ⎜− sin θ12 ⎜ ⎝ 0 ⎛ cos θ13 ⎜ ⎜ 0 ⎜ ⎝− sin θ13

(A.2c)

0 0 ⎞ sin θ12 0 ⎞⎛ 1 ⎟ ⎟⎜ θ θ23 ⎟ 0 cos sin 23 cos θ12 0 ⎟⎜ ⎟⎜ ⎟ 0 1 ⎠⎝ 0 − sin θ23 cos θ23 ⎠ ⎞ 0 sin θ13 ⎟ 1 0 ⎟ ⎟ 0 cos θ13 ⎠

⎛ cos θ12cos θ13 sin θ12cos θ23 cos θ12 sin θ13cos θ23 ⎞ ⎜ ⎟ ⎜ − cos θ12 sin θ13sin θ23 +sin θ12 sin θ23 ⎟ ⎜ ⎟ −sin θ12 sin θ13cos θ23 ⎟ cos θ12cos θ23 = ⎜− sin θ12 cos θ13 ⎜ ⎟ + sin θ12 sin θ13sin θ23 + cos θ12 sin θ23 ⎟ ⎜ ⎜ ⎟ − cos θ13sin θ23 cos θ13cos θ23 ⎝ − sin θ13 ⎠

⎛ cos θ12cos θ13 sin θ12cos θ23 cos θ12sin θ13 ⎞ ⎜ ⎟ ⎜ − sin θ12 sin θ13sin θ23 +sin θ12cos θ13sin θ23⎟ ⎜ ⎟ −sin θ12 sin θ13 ⎟ cos θ12cos θ23 = ⎜ − sin θ12cos θ13 ⎜ ⎟ + cos θ12cos θ13sin θ23⎟ ⎜− cos θ12 sin θ13sin θ23 ⎜ ⎟ − sin θ23 cos θ13cos θ23 ⎠ ⎝ − sin θ13cos θ23

(A.3)

The diabatic matrix elements as function of the ADT angles are as follows:

(A.1) 3492



Article

S.I: (x0 , y0 ) → (x0 + Δx , y0 ) → (x0 + Δx , y0 + Δy)

W11 = u1 cos2 θ12 cos2 θ13 2

+ u3(cos θ12 cos θ23 sin θ13 + sin θ12 sin θ23)

S.II: (x0 , y0 ) → (x0 , y0 + Δy) → (x0 + Δx , y0 + Δy)

2

+ u 2(cos θ23 sin θ12 − cos θ12 sin θ13 sin θ23)

The corresponding ADT matrices for both the paths and steps are given by

(A.4a)

W22 = u1 cos 2 θ13 sin 2 θ12

S.I: A(x0 + Δx , y0 ) = exp(−

+ u3( −cos θ23 sin θ12 sin θ13 + cos θ12 sin θ23)2

τx(x , y0 )dx)A(x0 , y0 )

0

A(x0 + Δx , y0 + Δy)

(A.4b)

W33 = u3 cos2 θ13 cos2 θ23 + u1 sin 2 θ13

= exp(−

y0 +Δy

∫y

τy(x0 + Δx , y)dy)A(x0 + Δx , y0 )

0

2

+ u 2 cos θ13 sin θ23

(A.4c)

(B.3b)

W12 = W21 = −u1 cos θ12 cos2 θ13 sin θ12

S.II: A(x0 , y0 + Δy)

+ u3[(− cos θ23 sin θ12 sin θ13 + cos θ12 sin θ23)

= exp(−

× (cos θ12 cos θ23 sin θ13 + sin θ12 sin θ23)] × (cos θ12 cos θ23 + sin θ12 sin θ13 sin θ23)]

y0 +Δy

∫y

τy(x0 , y)dy)A(x0 , y0 )

(B.4a)

0

+ u 2[(cos θ23 sin θ12 − cos θ12 sin θ13 sin θ23)

A(x0 + Δx , y0 + Δy) (A.4d)

= exp(−

∫x

W13 = W31

x0 +Δx 0

τx(x , y0 + Δy)dx)A(x0 , y0 + Δy) (B.4b)

When eqs B.3a and B.4a are substituted in eqs B.3b and B.4b, respectively, the corresponding ADT matrices for two different chosen paths will be

= cos θ13[−(u 2 − u3)cos θ23 sin θ12 sin θ23 2

2

+ cos θ12( −u1 + u3 cos θ23 + u 2 sin θ23)sin θ13] (A.4e)

S.I: A1 ≡A(x0 + Δx , y0 + Δy)

W23 = W32

= exp(−

= cos θ13[−(u 2 − u3)cos θ12 cos θ23 sin θ23

−

(A.4f)

■

∂ A(x , y) + τyA(x , y) = 0 ∂y

(B.1b)

τx(x , y0 )dx)A(x0 , y0 )

0

(B.5a)

∫y

∫x

x0 +Δx 0

y0 +Δy

τx(x , y0 + Δy)dx

τy(x0 , y)dy)A(x0 , y0 )

(B.5b)

Since the NAC matrix is skew-symmetric and its dagger undergoes sign change, we have A1† = A†(x0 , y0 ) exp(

∫y

y0 +Δy

τy(x0 + Δx , y)dy

0

+

s 0

τy(x0 + Δx , y)dy

0

as τ(s′|Γ)ds′)A(s0)

x0 +Δx

= exp(−

The solution for any of these scalar equations can be written

∫s

∫x

− (B.1a)

y0 +Δy

S.II: A 2 ≡A(x0 + Δx , y0 + Δy)

APPENDIX B: PATH-DEPENDENT ADIABATIC TO DIABATIC TRANSFORMATION MATRICES We consider the following adiabatic-to-diabatic transformation (ADT) equation characterized by the coordinates x and y: ∂ A(x , y) + τxA(x , y) = 0 ∂x

∫y

0

+ sin θ12(u1 − u3 cos2 θ23 − u 2 sin 2 θ23)sin θ13]

A(s|Γ) = ς exp( −

x0 +Δx

(B.3a)

+ u 2(cos θ12 cos θ13 + sin θ12 sin θ13 sin θ23)2

2

∫x

(B.2)

∫x

x0 +Δx 0

τx(x , y0 )dx)

(B.6)

and therefore,

where s ∈ x or y. A(s0)(= I) and A(s|Γ) are the orthogonal transformation matrices at s0 and s, respectively. The ordering operator ς dictates the order of integration along an assign contour, where the grid of L points (P0, P1,..., Pn,..., PL) define the contour, Γ. We carry out two consecutive steps of integration in two different paths Γ1 and Γ2 to reach a specific point and find the corresponding ADT matrices, A1(x, y) and A2(x, y) by employing eq B.2. The two different paths on the xy plane are chosen as follows:

x0 +Δx

A1†A 2 = A†(x0 , y0 ) exp( −

∫x

− τx(x , y0 )]dx +

∫y

0

y0 +Δy

[τx(x , y0 + Δy) [τy(x0 + Δx , y)

0

− τy(x0 , y)]dy)A(x0 , y0 )

(B.7a)

Considering a midpoint (x,̅ y)̅ in the planar interval as x0 ≤ x̅ ≤ x0 + Δx and y0 ≤ y ̅ ≤ y0 + Δy and assuming Δx,Δy are small 3493



Article

enough, we use the definition of numerical differentiation as well as integration around the midpoint (x,̅ y)̅ and safely write eq B.7a as: A1†A 2

div τ 12 ⃗ = 2 sin θ12 cos θ12 cos2 θ23(∇⃗θ13·∇⃗θ13) − 2 sin θ12 cos θ12(∇⃗θ23·∇⃗θ23) − 3 cos2 θ12 cos θ23 (∇⃗θ13·∇⃗θ23) + sin 2 θ12 cos θ23(∇⃗θ13·∇⃗θ23)

⎛ x0 +Δx ∂τ (x , y ) x ̅ dx =A (x0 , y0 ) exp⎜⎜ −Δy x0 ∂y ⎝ y0 +Δy ∂τy(x , y) ⎞ ̅ dy⎟⎟A(x0 , y0 ) + Δx y0 ∂x ⎠

∫

†

− sin θ23∇2 θ13 − ∇2 θ12

∫

div τ ⃗ 23 = 2 sin θ12 sin θ23 cos θ23(∇⃗θ13·∇⃗θ13) + 3 cos θ12 cos θ23(∇⃗θ12·∇⃗θ23) + 3 sin θ12(∇⃗θ12·∇⃗θ23)

⎛ ⎡ ∂τ (x , y ) = A†(x0 , y0 ) exp⎜⎜ −ΔxΔy⎢ x ̅ ̅ ⎢⎣ ∂y ⎝ ⎞ ∂τy(x ̅ , y ̅ ⎤ ⎥⎟A(x0 , y0 ) − ∂x ⎥⎦⎟⎠

+ sin θ12 sin θ23(∇⃗θ13·∇⃗θ23) + sin θ12 cos θ23∇⃗θ12 − cos θ12∇2 θ23

+ 3 sin θ12 cos θ23(∇⃗θ12·∇⃗θ13) − 3 cos θ12(∇⃗θ12·∇⃗θ23) 2

which in turn takes the following form:

− cos θ12 sin θ23(∇⃗θ13·∇⃗θ23) − cos θ12 cos θ23∇⃗ θ13

A1†A 2 = A†(x0 , y0 ) exp( −ΔxΔy(curl τ )xy )A(x0 , y0 )

■

(B.7c)

Since the argument of the exponential in eq B.7c is not a null matrix, the product matrix (B = A†1A2) does not appear as unit one, and thereby, the ADT matrices, A1 and A2, are not the same at the point (x0 + Δx; y0 + Δy) for two different chosen paths of integration, S.I and S.II. If the curl of the NAC matrix is zero (i.e., curl τ = 0), both the chosen path (even any path) will produce the same ADT matrices. On the contrary, the product matrix (B = A†1A2) is an orthogonal one as shown below: B† = A†(x0 , y0 ) exp(ΔxΔy(curl τ )xy )A(x0 , y0 )

■

AUTHOR INFORMATION

Corresponding Author

*(S.A.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

exp( −ΔxΔy(curl τ )xy )A(x0 , y0 )

■

(B.9)

=I

ACKNOWLEDGMENTS S.A. acknowledges Professor Koushik Roy, Department of Theoretical Physics, IACS, for stimulating discussion on the proof of the orthogonality of the path-dependent ADT matrices. S.M. acknowledges INSPIRE Programme (IF110026) of DST, India, for a research fellowship, and S.A. acknowledges DST, India, through project no. SR/S1/PC-13/ 2008 for the research funding.

†

where curl τ⃗ is a skew-symmetric matrix, (curl τ⃗) = −curl τ⃗. Therefore, the ADT matrices, A1 and A2, are related through an orthogonal transformation matrix, B.

■

APPENDIX C: CURL−DIVERGENCE EQUATIONS The explicit form of the curl equation14 in terms of ADT angles for each NACT is obtained by using eqs 10 and 28 as below:

■

12 curl τpq = C12 = Z12 = − cos θ23(∇q θ23∇p θ13 − ∇p θ23∇q θ13)

(C.1a)

REFERENCES

(1) Born, M.; Oppenheimer, J. R. Ann. Phys. 1927, 84, 457. (2) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: Oxford, U.K., 1954. (3) Baer, M.; Niedner-Schatteburg, G.; Toennies, J. J. Chem. Phys. 1989, 91, 4169. (4) Baer, M.; Liao, C. L.; Flesch, R. X. G. D.; Nourbaksh, S.; Ng, C. Y. J. Chem. Phys. 1990, 93, 4845. (5) Aguilon, F.; Sizun, M.; Sidis, V.; Billing, G. D.; Markovic, N. J. Chem. Phys. 1996, 104, 4530. (6) Last, I.; Gilibert, M.; Baer, M. J. Chem. Phys. 1997, 107, 1451. (7) Grimbert, D.; Sidis, V.; Cobut, V. J. Chem. Phys. 1998, 108, 6331. (8) Baer, R.; Charutz, D.; Kosloff, R.; Baer, M. J. Chem. Phys. 1996, 105, 9141. (9) Charutz, D. M.; Baer, R.; Baer, M. Chem. Phys. Lett. 1997, 265, 629. (10) Adhikari, S.; Billing, G. D. J. Chem. Phys. 1999, 111, 40. (11) Varandas, A. J. C.; Xu, Z. R. J. Chem. Phys. 2000, 112, 2121.

= C23 = Z 23 = cos θ12 cos θ23(∇q θ12∇p θ13

− ∇p θ12∇q θ13) − sin θ12 sin θ23(∇q θ23∇p θ13 − ∇p θ23∇q θ13) + sin θ12(∇q θ12∇p θ23 − ∇p θ12∇q θ23) (C.1b) 13 curl τpq = C13 = Z13 = sin θ12 cos θ23(∇q θ12∇p θ13

− ∇p θ12∇q θ13) + cos θ12 sin θ23(∇q θ23∇p θ13 − ∇p θ23∇q θ13) − cos θ12(∇q θ12∇p θ23 − ∇p θ12∇q θ23) (C.1c)

where the divergences of τ⃗ (eq 10) are given by:

ASSOCIATED CONTENT

Property of orthogonality between the two ADT matrices due to the choice of two different paths and the details of nuclear dynamics (e.g., the profiles of diabatic population dynamics and the contour plots of diabatic probability densities). This material is available free of charge via the Internet at http:// pubs.acs.org.

†

ij

(C.2c)

* Supporting Information

B B = A (x0 , y0 ) exp(ΔxΔy(curl τ )xy )A(x0 , y0 )A (x0 , y0 )

curl

− sin θ12∇2 θ23

S

(B.8)

†

23 τpq

(C.2b)

div τ 13 ⃗ = 2 sin θ12 sin θ23 cos θ23(∇⃗θ12·∇⃗θ12) (B.7b)

†

(C.2a)

17

3494



Article

(12) Baer, M.; Lin, S. H.; Alijah, A.; Adhikari, S.; Billing, G. D. Phys. Rev. A 2000, 62, 32506. (13) Adhikari, S.; Billing, G. D.; Alijah, A.; Lin, S. H.; Baer, M. Phys. Rev. A 2000, 62, 32507. (14) Sarkar, B.; Adhikari, S. J. Chem. Phys. 2006, 124, 074101. (15) Sarkar, B.; Adhikari, S. Indian J. Phys. 2007, 81 (9), 925. (16) Sarkar, B.; Adhikari, S. J. Phys. Chem. A 2008, 112, 9868. (17) Sarkar, B.; Adhikari, S. Int. J. Quantum Chem. 2009, 109, 650. (18) Paul, A. K.; Sardar, S.; Sarkar, B.; Adhikari, S. J. Chem. Phys. 2009, 131, 124312. (19) Paul, A. K.; Sarkar, B.; Adhikari, S. In Recent Advances in Spectroscopy; Chaudhuri, R. K., Mekkaden, M. V., Raveendran, A. V., Narayanan, A. S., Eds.; Astrophysics and Space Science Proceedings; Springer-Verlag: Berlin, Germany, 2010; Chapter 7, p 63. (20) Higgins, H. C. L. Adv. Spectrosc. 1961, 2, 429. (21) Higgins, H. C. L. Proc. R. Soc. London, Ser. A 1975, 344, 147. (22) Manolopoulos, D. E.; Child, M. S. Phys. Rev. Lett. 1999, 82, 2223. (23) Child, M. S. Adv. Chem. Phys. 2002, 124, 1. (24) Varandas, A. J. C.; Tennyson, J.; Murrel, J. N. Chem. Phys. Lett. 1979, 61, 431. (25) Herzberg, G.; Higgins, H. C. L. Discuss. Faraday Soc. 1963, 35, 77. (26) Hobey, W. D.; Mclachlan, A. D. J. Chem. Phys. 1960, 33, 1695. (27) Smith, F. T. Phys. Rev. 1969, 179, 111. (28) Baer, M. Chem. Phys. Lett. 1975, 35, 112. (29) Baer, M. Phys. Rep. 2002, 358, 75. (30) Baer, M. Beyond Born−Oppenheimer: Conical Intersections and Electronic Nonadiabatic Coupling Terms; Wiley-Interscience: New York, 2006. (31) Mead, C. A.; Truhlar, D. G. J. Chem. Phys. 1982, 77, 6090. (32) Kendrick, B. K.; Mead, C. A.; Truhlar, D. G. Chem. Phys. 2002, 277, 31. (33) Matsunaga, N.; Yarkony, D. R. Mol. Phys. 1998, 93, 79. (34) Baer, M.; Englman, R. Mol. Phys. 1992, 75, 293. (35) Alijah, A.; Baer, M. J. Phys. Chem. A 2000, 104, 389. (36) Yarkony, D. R. J. Chem. Phys. 1996, 105, 10456. (37) Sadygov, R. G.; Yarkony, D. R. J. Chem. Phys. 1998, 109, 20. (38) Vértesi, T.; Vibók, A.; Halász, G. J.; Baer, M. J. Chem. Phys. 2004, 121, 4000. (39) Abrol, R.; Kuppermann, A. J. Chem. Phys. 2002, 116, 1035. (40) Hellmann, H. Einfuhrang in die Quantenchemie; Franz Deutiche: Leipzig, Germany, 1937. (41) Feynman, R. Phys. Rev. 1939, 56, 340. (42) Epstein, S. T. Am. J. Phys. 1954, 22, 613. (43) Sadygov, R. G.; Yarkony, D. R. J. Chem. Phys. 1999, 110, 3639. (44) Koizumi, H.; Bersuker, I. B. Phys. Rev. Lett. 1999, 83, 3009. (45) Cocchini, F.; Upton, T. H.; Andreoni, W. J. Chem. Phys. 1988, 88, 6068. (46) Paul, A. K.; Ray, S.; Mukhopadhyay, D.; Adhikari, S. Chem. Phys. Lett. 2011, 508, 300. (47) Paul, A. K.; Ray, S.; Mukhopadhyay, D.; Adhikari, S. J. Chem. Phys. 2011, 135, 034107. (48) Paul, A. K.; Ray, S.; Adhikari, S. Vibronic Interactions and the Jahn−Teller Effect: Theory and Applications; Progress in Theoretical Chemistry and Physics 23; Springer: New York, 2012; Chapter 15, p 281. (49) Delacrétaz, G.; Wöste, L. Surf. Sci. 1985, 156, 770. (50) Delacrétaz, G.; Grant, E. R.; Whetten, R. L.; Wöste, L.; Zwanziger, J. W. Phys. Rev. Lett. 1986, 56, 2598. (51) Rutz, S.; Kobe, K.; Kühling, H.; Schreiber, E.; Wöste, L. Z. Phys. D 1993, 26, 276. (52) Ernst, W. E.; Rakowsky, S. Phys. Rev. Lett. 1995, 74, 58. (53) Ohashi, N.; Tsuura, M.; Hougen, J. T.; Ernst, W. E.; Rakowsky, S. J. Mol. Spectrosc. 1997, 184, 22. (54) Ernst, W. E.; Golonzka, O. Phys. Scr. 2004, T112, 27. (55) Higgins, H. C. L. Mol. Phys. 1963, 6, 445. (56) Al-Jabour, S.; Baer, M.; Deeb, O.; Leibscher, M.; Manz, J.; Xu, X.; Zilberg, S. J. Phys. Chem. A 2010, 114, 2991.

(57) Top, Z. H.; Baer, M. J. Chem. Phys. 1977, 66, 1363. (58) Bunker, P. R.; Jensen, P. Molecular Symmetry and Spectroscopy, 2nd ed.; National Research Council of Canada: Ontario, Canada, 2006. (59) Kettle, S. F. A. Symmetry and Structure: Readable Group Theory for Chemists, 3rd ed.; J. Wiley and Sons Inc: West Sussex, U.K., 2007. (60) Vértesi, T.; Vibók, A.; Halász, G. J.; Baer, M. J. Chem. Phys. 2004, 120, 2565. (61) Werner, H.-J.; et al. MOLPRO, version 2010.1, a package of ab initio programs, 2010; see http://www.molpro.net. (62) Vibók, A.; Halász, G. J.; Mebel, A. M.; Hu, S.; Baer, M. Int. J. Quantum Chem. 2004, 99, 594. (63) Abrol, R.; Shaw, A.; Kupperman, A.; Yarkony, D. R. J. Chem. Phys. 2001, 115, 4640. (64) Baer, M. Chem. Phys. Lett. 2000, 329, 450. (65) Meiswinkel, R.; Kóppel, H. Chem. Phys. 1990, 144, 117. (66) de Vivie-Riedle, R.; Gaus, J.; Bonačić-Koutecký, V.; Manz, J.; Reischl-Lenz, B.; Saalfrank, P. Chem. Phys. 1997, 223, 1.

3495


Construction of Diabatic Hamiltonian Matrix from ab Initio Calculated

Recommend Documents